diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 6f290b6a..60029272 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -24,7 +24,7 @@ jobs: steps: - name: Check out repository uses: actions/checkout@v2 - + - name: Build environment uses: conda-incubator/setup-miniconda@v2 with: @@ -32,7 +32,7 @@ jobs: miniforge-variant: Mambaforge miniforge-version: 4.9.2-4 use-mamba: true - + - name: Install package shell: bash -l {0} run: pip install -e . @@ -40,7 +40,7 @@ jobs: - name: Run pytest with coverage report shell: bash -l {0} run: python -m pytest -rs -v --cov=./ --cov-report=xml - + - name: Upload coverage to Codecov uses: codecov/codecov-action@v3 with: diff --git a/.gitignore b/.gitignore index b6e47617..b11f06ba 100644 --- a/.gitignore +++ b/.gitignore @@ -127,3 +127,6 @@ dmypy.json # Pyre type checker .pyre/ + +# IDE files +.idea/ diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 00000000..11025bc5 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,47 @@ +ci: + autofix_commit_msg: | + [pre-commit.ci] auto fixes from pre-commit.com hooks + + for more information, see https://pre-commit.ci + autofix_prs: true + autoupdate_commit_msg: '[pre-commit.ci] pre-commit autoupdate' + autoupdate_schedule: weekly + skip: [] + submodules: false + +repos: +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.4.0 + hooks: + - id: check-yaml + - id: end-of-file-fixer + - id: trailing-whitespace + exclude: 'hoomd_polymers/tests/assets/.* | hoomd_polymers/assets/.*' +- repo: https://github.com/psf/black + rev: 23.7.0 + hooks: + - id: black + args: [--line-length=80] +- repo: https://github.com/pycqa/isort + rev: 5.12.0 + hooks: + - id: isort + name: isort (python) + args: + [--profile=black, --line-length=80] + exclude: 'cmeuitls/tests/assets/.* ' + +- repo: https://github.com/pycqa/flake8 + rev: 6.1.0 + hooks: + - id: flake8 + args: + - --max-line-length=80 + exclude: '__init__.py' + +#- repo: https://github.com/pycqa/pydocstyle +# rev: '6.3.0' +# hooks: +# - id: pydocstyle +# exclude: ^(cmeuitls/tests/|setup.py) +# args: [--convention=numpy] diff --git a/README.md b/README.md index 66017fd5..174fb692 100644 --- a/README.md +++ b/README.md @@ -2,25 +2,25 @@ [![codecov](https://codecov.io/gh/chrisjonesBSU/hoomd-polymers/branch/main/graph/badge.svg?token=86LY9WHSH6)](https://codecov.io/gh/chrisjonesBSU/hoomd-polymers) ## About -Wrapper for [MoSDeF](https://github.com/mosdef-hub) packages and [Hoomd-Blue](https://github.com/glotzerlab/hoomd-blue) with +Wrapper for [MoSDeF](https://github.com/mosdef-hub) packages and [Hoomd-Blue](https://github.com/glotzerlab/hoomd-blue) with a focus on simulating soft matter systems. ## Installation -### 1. Clone this repository: ### +### 1. Clone this repository: ### ``` -git clone git@github.com:chrisjonesbsu/hoomd-polymers.git -cd hoomd-polymers +git clone git@github.com:chrisjonesbsu/hoomd-polymers.git +cd hoomd-polymers ``` -### 2. Set up and activate environment: ### +### 2. Set up and activate environment: ### #### a. Using HOOMD-blue from conda: ``` -conda env create -f environment-cpu.yml -conda activate hoomd-polymers +conda env create -f environment-cpu.yml +conda activate hoomd-polymers python -m pip install . -``` +``` ## Basic Usage #### Using the built in molecules, systems and forcefields: diff --git a/codecov.yml b/codecov.yml index 8cf18d8d..79f944f5 100644 --- a/codecov.yml +++ b/codecov.yml @@ -8,7 +8,7 @@ comment: - "main" ignore: - "hoomd-polymers/tests" - - "hoomd-polymers/library/*" - - "hoomd-polymers/logs" + - "hoomd-polymers/assets" - "hoomd-polymers/__version__.py" + - "tutorials" - "setup.py" diff --git a/environment-cpu.yml b/environment-cpu.yml index a57340a1..4e2ced7a 100644 --- a/environment-cpu.yml +++ b/environment-cpu.yml @@ -1,19 +1,26 @@ -name: hoomd-polymers +name: hoomd-polymers channels: - conda-forge dependencies: + - boltons - foyer - freud - - gsd - - hoomd>3.7=*cpu* + # - gmso>=0.11.0 + - gsd<3.0 + - hoomd=3.11=*cpu* - mbuild - numpy + - openbabel - pip + - py3Dmol + - pydantic - pytest - pytest-cov - python=3.9 - - signac=1.7.0 - - signac-flow=0.18.1 + - symengine + - python-symengine - unyt - pip: - git+https://github.com/cmelab/cmeutils.git@master + - git+https://github.com/cmelab/grits.git@main + - git+https://github.com/mosdef-hub/gmso@main diff --git a/environment-dev.yml b/environment-dev.yml new file mode 100644 index 00000000..c375785c --- /dev/null +++ b/environment-dev.yml @@ -0,0 +1,23 @@ +name: hoomd-polymers-dev +channels: + - conda-forge +dependencies: + - foyer + - freud + #- gmso >= 0.11.0 + - gsd<3.0 + - hoomd=3.11=*gpu* + - mbuild + - numpy + - openbabel + - pip + - py3Dmol + - pytest + - pytest-cov + - python=3.9 + - unyt + - pre-commit + - pip: + - git+https://github.com/cmelab/cmeutils.git@master + - git+https://github.com/cmelab/grits.git@main + - git+https://github.com/mosdef-hub/gmso@main diff --git a/environment-gpu.yml b/environment-gpu.yml index e757992e..0cf8a2ca 100644 --- a/environment-gpu.yml +++ b/environment-gpu.yml @@ -1,19 +1,26 @@ -name: hoomd-polymers +name: hoomd-polymers channels: - conda-forge dependencies: + - boltons - foyer - freud - - gsd - - hoomd>3.0=*gpu* + # - gmso>=0.11.0 + - gsd<3.0 + - hoomd=3.11=*gpu* - mbuild - numpy + - openbabel - pip + - py3Dmol + - pydantic - pytest - pytest-cov - python=3.9 - - signac=1.7.0 - - signac-flow=0.18.1 + - symengine + - python-symengine - unyt - pip: - git+https://github.com/cmelab/cmeutils.git@master + - git+https://github.com/cmelab/grits.git@main + - git+https://github.com/mosdef-hub/gmso@main diff --git a/environment-nohoomd.yml b/environment-nohoomd.yml deleted file mode 100644 index c66c25d8..00000000 --- a/environment-nohoomd.yml +++ /dev/null @@ -1,18 +0,0 @@ -name: hoomd-polymers -channels: - - conda-forge -dependencies: - - foyer - - freud - - gsd - - mbuild - - numpy - - pip=22.2 - - pytest - - pytest-cov - - python=3.9 - - signac=1.7.0 - - signac-flow=0.18.1 - - unyt - - pip: - - git+https://github.com/cmelab/cmeutils.git@master diff --git a/hoomd_polymers/__init__.py b/hoomd_polymers/__init__.py index 4483c6a2..25815e39 100644 --- a/hoomd_polymers/__init__.py +++ b/hoomd_polymers/__init__.py @@ -1,2 +1,9 @@ -from hoomd_polymers.sim import simulation -from hoomd_polymers.modules import welding +from .base import ( + CoPolymer, + Lattice, + Molecule, + Pack, + Polymer, + Simulation, + System, +) diff --git a/hoomd_polymers/assets/__init__.py b/hoomd_polymers/assets/__init__.py new file mode 100644 index 00000000..00f65054 --- /dev/null +++ b/hoomd_polymers/assets/__init__.py @@ -0,0 +1,2 @@ +from .forcefields import FF_DIR +from .molecule_files import MON_DIR diff --git a/hoomd_polymers/library/forcefields/__init__.py b/hoomd_polymers/assets/forcefields/__init__.py similarity index 98% rename from hoomd_polymers/library/forcefields/__init__.py rename to hoomd_polymers/assets/forcefields/__init__.py index 37798e21..7c6deae0 100644 --- a/hoomd_polymers/library/forcefields/__init__.py +++ b/hoomd_polymers/assets/forcefields/__init__.py @@ -1,4 +1,3 @@ import os - FF_DIR = os.path.abspath(os.path.dirname(__file__)) diff --git a/hoomd_polymers/assets/forcefields/benzene_opls.xml b/hoomd_polymers/assets/forcefields/benzene_opls.xml new file mode 100644 index 00000000..822b6423 --- /dev/null +++ b/hoomd_polymers/assets/forcefields/benzene_opls.xml @@ -0,0 +1,23 @@ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/hoomd_polymers/assets/forcefields/dimethylether_opls.xml b/hoomd_polymers/assets/forcefields/dimethylether_opls.xml new file mode 100644 index 00000000..6909424b --- /dev/null +++ b/hoomd_polymers/assets/forcefields/dimethylether_opls.xml @@ -0,0 +1,25 @@ + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/hoomd_polymers/library/forcefields/gaff.xml b/hoomd_polymers/assets/forcefields/gaff.xml similarity index 100% rename from hoomd_polymers/library/forcefields/gaff.xml rename to hoomd_polymers/assets/forcefields/gaff.xml diff --git a/hoomd_polymers/assets/forcefields/oplsaa.xml b/hoomd_polymers/assets/forcefields/oplsaa.xml new file mode 100644 index 00000000..b05d8527 --- /dev/null +++ b/hoomd_polymers/assets/forcefields/oplsaa.xml @@ -0,0 +1,3948 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/hoomd_polymers/library/forcefields/pps_opls.xml b/hoomd_polymers/assets/forcefields/pps_opls.xml similarity index 100% rename from hoomd_polymers/library/forcefields/pps_opls.xml rename to hoomd_polymers/assets/forcefields/pps_opls.xml diff --git a/hoomd_polymers/assets/molecule_files/IPH.mol2 b/hoomd_polymers/assets/molecule_files/IPH.mol2 new file mode 100644 index 00000000..a861d588 --- /dev/null +++ b/hoomd_polymers/assets/molecule_files/IPH.mol2 @@ -0,0 +1,34 @@ +@MOLECULE +***** + 13 13 0 0 0 +SMALL +GASTEIGER + +@ATOM + 1 C 0.0000 -0.0200 0.9120 C.ar 1 UNL1 0.1169 + 2 C 1.2000 -0.0130 0.2170 C.ar 1 UNL1 -0.0202 + 3 C 1.1980 0.0040 -1.1630 C.ar 1 UNL1 -0.0583 + 4 C 0.0000 0.0130 -1.8540 C.ar 1 UNL1 -0.0615 + 5 C -1.1980 0.0070 -1.1630 C.ar 1 UNL1 -0.0583 + 6 C -1.2000 -0.0160 0.2170 C.ar 1 UNL1 -0.0202 + 7 O 0.0000 -0.0370 2.2710 O.3 1 UNL1 -0.5068 + 8 H 2.1360 -0.0200 0.7560 H 1 UNL1 0.0654 + 9 H 2.1320 0.0100 -1.7050 H 1 UNL1 0.0619 + 10 H 0.0000 0.0270 -2.9340 H 1 UNL1 0.0618 + 11 H -2.1320 0.0150 -1.7050 H 1 UNL1 0.0619 + 12 H -2.1360 -0.0210 0.7550 H 1 UNL1 0.0654 + 13 H 0.0000 0.8840 2.5620 H 1 UNL1 0.2921 +@BOND + 1 1 2 ar + 2 1 6 ar + 3 1 7 1 + 4 2 3 ar + 5 2 8 1 + 6 3 4 ar + 7 3 9 1 + 8 4 5 ar + 9 4 10 1 + 10 5 6 ar + 11 5 11 1 + 12 6 12 1 + 13 7 13 1 diff --git a/hoomd_polymers/library/monomers/__init__.py b/hoomd_polymers/assets/molecule_files/__init__.py similarity index 98% rename from hoomd_polymers/library/monomers/__init__.py rename to hoomd_polymers/assets/molecule_files/__init__.py index a4526919..1d0610d3 100644 --- a/hoomd_polymers/library/monomers/__init__.py +++ b/hoomd_polymers/assets/molecule_files/__init__.py @@ -1,4 +1,3 @@ import os - MON_DIR = os.path.abspath(os.path.dirname(__file__)) diff --git a/hoomd_polymers/library/monomers/pekk_meta.mol2 b/hoomd_polymers/assets/molecule_files/pekk_meta.mol2 similarity index 89% rename from hoomd_polymers/library/monomers/pekk_meta.mol2 rename to hoomd_polymers/assets/molecule_files/pekk_meta.mol2 index d8adb428..eca32717 100644 --- a/hoomd_polymers/library/monomers/pekk_meta.mol2 +++ b/hoomd_polymers/assets/molecule_files/pekk_meta.mol2 @@ -6,43 +6,43 @@ NO_CHARGES @CRYSIN 9.8510 10.2351 20.1594 90.0000 90.0000 90.0000 1 1 @ATOM - 1 C 0.2404 -0.8324 0.8626 C 1 RES - 2 C -0.1273 -1.4670 -0.3286 C 1 RES - 3 C -0.4484 -0.7059 -1.4570 C 1 RES - 4 C -0.4020 0.6904 -1.3952 C 1 RES - 5 C -0.0345 1.3268 -0.2048 C 1 RES - 6 C 0.2876 0.5683 0.9271 C 1 RES - 7 O 0.6521 1.2133 2.1067 O 1 RES - 8 H 0.4877 -1.4277 1.7323 H 1 RES - 9 H -0.1633 -2.5476 -0.3770 H 1 RES - 10 H -0.6506 1.2790 -2.2687 H 1 RES - 11 H 0.0000 2.4079 -0.1616 H 1 RES - 12 C 0.9526 0.2064 5.7545 C 1 RES - 13 C 0.6242 0.9234 4.5994 C 1 RES - 14 C 1.0055 0.4425 3.3441 C 1 RES - 15 C 1.7160 -0.7562 3.2413 C 1 RES - 16 C 2.0468 -1.4769 4.3936 C 1 RES - 17 C 1.6664 -0.9995 5.6611 C 1 RES - 18 C 2.0046 -1.7428 6.8886 C 1 RES - 19 O 2.6439 -2.8272 6.8372 O 1 RES - 20 H 0.6496 0.5926 6.7202 H 1 RES - 21 H 0.0739 1.8521 4.6778 H 1 RES - 22 H 2.0108 -1.1277 2.2684 H 1 RES - 23 H 2.5976 -2.4037 4.2932 H 1 RES - 24 C 0.8032 -0.2154 10.7215 C 1 RES - 25 C -0.0000 0.0000 9.5970 C 1 RES - 26 C 0.3911 -0.4969 8.3508 C 1 RES - 27 C 1.5863 -1.2102 8.2270 C 1 RES - 28 C 2.3929 -1.4282 9.3487 C 1 RES - 29 C 2.0061 -0.9311 10.6066 C 1 RES - 30 C 2.8411 -1.1493 11.8020 C 1 RES - 31 O 3.9251 -1.7872 11.7308 O 1 RES - 32 H 0.4860 0.1764 11.6803 H 1 RES - 33 H -0.9259 0.5523 9.6917 H 1 RES - 34 H -0.2316 -0.3294 7.4817 H 1 RES - 35 H 3.3154 -1.9832 9.2321 H 1 RES - 36 H -0.7358 -1.2029 -2.3883 H 1 RES - 37 H 2.5223 -0.7544 12.7712 H 1 RES + 1 C 0.2404 -0.8324 0.8626 C 1 RES + 2 C -0.1273 -1.4670 -0.3286 C 1 RES + 3 C -0.4484 -0.7059 -1.4570 C 1 RES + 4 C -0.4020 0.6904 -1.3952 C 1 RES + 5 C -0.0345 1.3268 -0.2048 C 1 RES + 6 C 0.2876 0.5683 0.9271 C 1 RES + 7 O 0.6521 1.2133 2.1067 O 1 RES + 8 H 0.4877 -1.4277 1.7323 H 1 RES + 9 H -0.1633 -2.5476 -0.3770 H 1 RES + 10 H -0.6506 1.2790 -2.2687 H 1 RES + 11 H 0.0000 2.4079 -0.1616 H 1 RES + 12 C 0.9526 0.2064 5.7545 C 1 RES + 13 C 0.6242 0.9234 4.5994 C 1 RES + 14 C 1.0055 0.4425 3.3441 C 1 RES + 15 C 1.7160 -0.7562 3.2413 C 1 RES + 16 C 2.0468 -1.4769 4.3936 C 1 RES + 17 C 1.6664 -0.9995 5.6611 C 1 RES + 18 C 2.0046 -1.7428 6.8886 C 1 RES + 19 O 2.6439 -2.8272 6.8372 O 1 RES + 20 H 0.6496 0.5926 6.7202 H 1 RES + 21 H 0.0739 1.8521 4.6778 H 1 RES + 22 H 2.0108 -1.1277 2.2684 H 1 RES + 23 H 2.5976 -2.4037 4.2932 H 1 RES + 24 C 0.8032 -0.2154 10.7215 C 1 RES + 25 C -0.0000 0.0000 9.5970 C 1 RES + 26 C 0.3911 -0.4969 8.3508 C 1 RES + 27 C 1.5863 -1.2102 8.2270 C 1 RES + 28 C 2.3929 -1.4282 9.3487 C 1 RES + 29 C 2.0061 -0.9311 10.6066 C 1 RES + 30 C 2.8411 -1.1493 11.8020 C 1 RES + 31 O 3.9251 -1.7872 11.7308 O 1 RES + 32 H 0.4860 0.1764 11.6803 H 1 RES + 33 H -0.9259 0.5523 9.6917 H 1 RES + 34 H -0.2316 -0.3294 7.4817 H 1 RES + 35 H 3.3154 -1.9832 9.2321 H 1 RES + 36 H -0.7358 -1.2029 -2.3883 H 1 RES + 37 H 2.5223 -0.7544 12.7712 H 1 RES @BOND 1 2 1 1 2 6 1 1 diff --git a/hoomd_polymers/library/monomers/pekk_para.mol2 b/hoomd_polymers/assets/molecule_files/pekk_para.mol2 similarity index 89% rename from hoomd_polymers/library/monomers/pekk_para.mol2 rename to hoomd_polymers/assets/molecule_files/pekk_para.mol2 index 987c7237..5fba82cd 100644 --- a/hoomd_polymers/library/monomers/pekk_para.mol2 +++ b/hoomd_polymers/assets/molecule_files/pekk_para.mol2 @@ -6,43 +6,43 @@ NO_CHARGES @CRYSIN 11.2871 9.9340 19.9925 90.0000 90.0000 90.0000 1 1 @ATOM - 1 C -0.0622 -0.6814 1.0132 C 1 RES - 2 C 0.0330 -1.5040 -0.1142 C 1 RES - 3 C 0.1160 -0.9367 -1.3898 C 1 RES - 4 C 0.1040 0.4537 -1.5389 C 1 RES - 5 C 0.0089 1.2780 -0.4126 C 1 RES - 6 C -0.0744 0.7138 0.8661 C 1 RES - 7 O -0.1687 1.5448 1.9801 O 1 RES - 8 H -0.1261 -1.1273 1.9976 H 1 RES - 9 H 0.0423 -2.5802 0.0006 H 1 RES - 10 H 0.1683 0.8923 -2.5262 H 1 RES - 11 H -0.0000 2.3538 -0.5328 H 1 RES - 12 C 0.5176 0.9298 5.6627 C 1 RES - 13 C 0.6009 1.4496 4.3669 C 1 RES - 14 C -0.2601 0.9857 3.3690 C 1 RES - 15 C -1.2063 0.0010 3.6649 C 1 RES - 16 C -1.2936 -0.5219 4.9593 C 1 RES - 17 C -0.4302 -0.0601 5.9695 C 1 RES - 18 C -0.5050 -0.5974 7.3404 C 1 RES - 19 O -1.3443 -1.4856 7.6462 O 1 RES - 20 H 1.1922 1.3000 6.4250 H 1 RES - 21 H 1.3335 2.2124 4.1366 H 1 RES - 22 H -1.8728 -0.3575 2.8911 H 1 RES - 23 H -2.0335 -1.2843 5.1688 H 1 RES - 24 C 0.8684 0.4761 10.7126 C 1 RES - 25 C -0.0000 0.0000 9.7251 C 1 RES - 26 C 0.4282 -0.0864 8.3978 C 1 RES - 27 C 1.7258 0.3031 8.0556 C 1 RES - 28 C 2.5979 0.7799 9.0396 C 1 RES - 29 C 2.1743 0.8700 10.3781 C 1 RES - 30 C 3.0770 1.3697 11.4313 C 1 RES - 31 O 4.2537 1.7303 11.1625 O 1 RES - 32 H 0.5212 0.5371 11.7369 H 1 RES - 33 H -1.0052 -0.3019 9.9891 H 1 RES - 34 H 2.0563 0.2357 7.0271 H 1 RES - 35 H 3.5994 1.0769 8.7544 H 1 RES - 36 H 0.1903 -1.5807 -2.2710 H 1 RES - 37 H 2.7281 1.4318 12.4663 H 1 RES + 1 C -0.0622 -0.6814 1.0132 C 1 RES + 2 C 0.0330 -1.5040 -0.1142 C 1 RES + 3 C 0.1160 -0.9367 -1.3898 C 1 RES + 4 C 0.1040 0.4537 -1.5389 C 1 RES + 5 C 0.0089 1.2780 -0.4126 C 1 RES + 6 C -0.0744 0.7138 0.8661 C 1 RES + 7 O -0.1687 1.5448 1.9801 O 1 RES + 8 H -0.1261 -1.1273 1.9976 H 1 RES + 9 H 0.0423 -2.5802 0.0006 H 1 RES + 10 H 0.1683 0.8923 -2.5262 H 1 RES + 11 H -0.0000 2.3538 -0.5328 H 1 RES + 12 C 0.5176 0.9298 5.6627 C 1 RES + 13 C 0.6009 1.4496 4.3669 C 1 RES + 14 C -0.2601 0.9857 3.3690 C 1 RES + 15 C -1.2063 0.0010 3.6649 C 1 RES + 16 C -1.2936 -0.5219 4.9593 C 1 RES + 17 C -0.4302 -0.0601 5.9695 C 1 RES + 18 C -0.5050 -0.5974 7.3404 C 1 RES + 19 O -1.3443 -1.4856 7.6462 O 1 RES + 20 H 1.1922 1.3000 6.4250 H 1 RES + 21 H 1.3335 2.2124 4.1366 H 1 RES + 22 H -1.8728 -0.3575 2.8911 H 1 RES + 23 H -2.0335 -1.2843 5.1688 H 1 RES + 24 C 0.8684 0.4761 10.7126 C 1 RES + 25 C -0.0000 0.0000 9.7251 C 1 RES + 26 C 0.4282 -0.0864 8.3978 C 1 RES + 27 C 1.7258 0.3031 8.0556 C 1 RES + 28 C 2.5979 0.7799 9.0396 C 1 RES + 29 C 2.1743 0.8700 10.3781 C 1 RES + 30 C 3.0770 1.3697 11.4313 C 1 RES + 31 O 4.2537 1.7303 11.1625 O 1 RES + 32 H 0.5212 0.5371 11.7369 H 1 RES + 33 H -1.0052 -0.3019 9.9891 H 1 RES + 34 H 2.0563 0.2357 7.0271 H 1 RES + 35 H 3.5994 1.0769 8.7544 H 1 RES + 36 H 0.1903 -1.5807 -2.2710 H 1 RES + 37 H 2.7281 1.4318 12.4663 H 1 RES @BOND 1 2 1 1 2 6 1 1 diff --git a/hoomd_polymers/base/__init__.py b/hoomd_polymers/base/__init__.py new file mode 100644 index 00000000..ee8a00ea --- /dev/null +++ b/hoomd_polymers/base/__init__.py @@ -0,0 +1,3 @@ +from .molecule import CoPolymer, Molecule, Polymer +from .simulation import Simulation +from .system import Lattice, Pack, System diff --git a/hoomd_polymers/base/molecule.py b/hoomd_polymers/base/molecule.py new file mode 100644 index 00000000..eb419f72 --- /dev/null +++ b/hoomd_polymers/base/molecule.py @@ -0,0 +1,417 @@ +import itertools +import os.path +import random +from typing import List + +import mbuild as mb +from gmso.core.topology import Topology +from gmso.external.convert_mbuild import from_mbuild, to_mbuild +from grits import CG_Compound +from mbuild.lib.recipes import Polymer as mbPolymer + +from hoomd_polymers.utils import check_return_iterable +from hoomd_polymers.utils.base_types import FF_Types +from hoomd_polymers.utils.exceptions import MoleculeLoadError +from hoomd_polymers.utils.ff_utils import ( + _validate_hoomd_ff, + apply_xml_ff, + find_xml_ff, +) + + +class Molecule: + def __init__( + self, num_mols, force_field=None, smiles=None, file=None, compound=None + ): + self.n_mols = num_mols + self.force_field = force_field + self.smiles = smiles + self.file = file + self.compound = compound + self._mapping = None + self._mb_molecule = self._load() + self._molecules = [] + self._cg_molecules = [] + self._generate() + self.gmso_molecule = self._convert_to_gmso(self._molecules[0]) + self._identify_topology_information(self.gmso_molecule) + if self.force_field: + self._validate_force_field() + + @property + def molecules(self): + """List of all instances of the molecule""" + if self._cg_molecules: + return self._cg_molecules + return self._molecules + + @property + def mapping(self): + """Dictionary of particle index to bead mapping""" + return self._mapping + + @mapping.setter + def mapping(self, mapping_array): + self._mapping = mapping_array + + @property + def n_particles(self): + n_particles = 0 + for molecule in self.molecules: + n_particles += molecule.n_particles + return n_particles + + @property + def n_bonds(self): + n_bonds = 0 + for molecule in self.molecules: + n_bonds += molecule.n_bonds + return n_bonds + + @property + def topology_information(self): + topology_information = dict() + topology_information["particle_types"] = self.particle_types + topology_information["hydrogen_types"] = self.hydrogen_types + topology_information["particle_charge"] = self.particle_charge + topology_information["particle_typeid"] = self.particle_typeid + topology_information["pair_types"] = self.pairs + topology_information["bond_types"] = self.bond_types + topology_information["angle_types"] = self.angle_types + topology_information["dihedral_types"] = self.dihedral_types + topology_information["improper_types"] = self.improper_types + return topology_information + + def coarse_grain(self, beads=None): + for comp in self.molecules: + cg_comp = CG_Compound(comp, beads=beads) + if cg_comp.mapping: + self._cg_molecules.append(cg_comp) + else: + raise ValueError( + "Unable to coarse grain the molecule. " + "Please check the bead types." + ) + if self._cg_molecules: + self.gmso_molecule = self._convert_to_gmso(self._cg_molecules[0]) + self._identify_topology_information(self.gmso_molecule) + + def _load(self): + if self.compound: + if isinstance(self.compound, mb.Compound): + return mb.clone(mb.clone(self.compound)) + if isinstance(self.compound, Topology): + return to_mbuild(self.compound) + else: + raise MoleculeLoadError( + msg=f"Unsupported compound type {type(self.compound)}. " + f"Supported compound types are: {str(mb.Compound)}" + ) + if self.file: + if isinstance(self.file, str) and os.path.isfile(self.file): + return mb.load(self.file) + else: + raise MoleculeLoadError( + msg=f"Unable to load the molecule from file {self.file}." + ) + + if self.smiles: + if isinstance(self.smiles, str): + return mb.load(self.smiles, smiles=True) + else: + raise MoleculeLoadError( + msg=f"Unable to load the molecule from smiles " + f"{self.smiles}." + ) + + def _generate(self): + for i in range(self.n_mols): + self._molecules.append(self._load()) + + def _convert_to_gmso(self, mb_molecule): + topology = from_mbuild(mb_molecule) + topology.identify_connections() + return topology + + def _identify_particle_information(self, gmso_molecule): + self.particle_types = [] + self.hydrogen_types = [] + self.particle_typeid = [] + self.particle_charge = [] + for site in gmso_molecule.sites: + p_name = getattr(site.atom_type, "name", None) or site.name + if p_name not in self.particle_types: + self.particle_types.append(p_name) + if ( + site.element + and site.element.atomic_number == 1 + and p_name not in self.hydrogen_types + ): + self.hydrogen_types.append(p_name) + self.particle_typeid.append(self.particle_types.index(p_name)) + self.particle_charge.append( + site.charge.to_value() if site.charge else 0 + ) + + def _identify_pairs(self, particle_types): + self.pairs = set( + itertools.combinations_with_replacement(particle_types, 2) + ) + + def _identify_bond_types(self, gmso_molecule): + self.bond_types = set() + for bond in gmso_molecule.bonds: + p1_name = ( + getattr(bond.connection_members[0].atom_type, "name", None) + or bond.connection_members[0].name + ) + p2_name = ( + getattr(bond.connection_members[1].atom_type, "name", None) + or bond.connection_members[1].name + ) + bond_connections = [p1_name, p2_name] + if not tuple(bond_connections[::-1]) in self.bond_types: + self.bond_types.add(tuple(bond_connections)) + + def _identify_angle_types(self, gmso_molecule): + self.angle_types = set() + for angle in gmso_molecule.angles: + p1_name = ( + getattr(angle.connection_members[0].atom_type, "name", None) + or angle.connection_members[0].name + ) + p2_name = ( + getattr(angle.connection_members[1].atom_type, "name", None) + or angle.connection_members[1].name + ) + p3_name = ( + getattr(angle.connection_members[2].atom_type, "name", None) + or angle.connection_members[2].name + ) + angle_connections = [p1_name, p2_name, p3_name] + if not tuple(angle_connections[::-1]) in self.angle_types: + self.angle_types.add(tuple(angle_connections)) + + def _identify_dihedral_types(self, gmso_molecule): + self.dihedral_types = set() + for dihedral in gmso_molecule.dihedrals: + p1_name = ( + getattr(dihedral.connection_members[0].atom_type, "name", None) + or dihedral.connection_members[0].name + ) + p2_name = ( + getattr(dihedral.connection_members[1].atom_type, "name", None) + or dihedral.connection_members[1].name + ) + p3_name = ( + getattr(dihedral.connection_members[2].atom_type, "name", None) + or dihedral.connection_members[2].name + ) + p4_name = ( + getattr(dihedral.connection_members[3].atom_type, "name", None) + or dihedral.connection_members[3].name + ) + dihedral_connections = [p1_name, p2_name, p3_name, p4_name] + if not tuple(dihedral_connections[::-1]) in self.dihedral_types: + self.dihedral_types.add(tuple(dihedral_connections)) + + def _identify_improper_types(self, gmso_molecule): + self.improper_types = set() + for improper in gmso_molecule.impropers: + p1_name = ( + getattr(improper.connection_members[0].atom_type, "name", None) + or improper.connection_members[0].name + ) + p2_name = ( + getattr(improper.connection_members[1].atom_type, "name", None) + or improper.connection_members[1].name + ) + p3_name = ( + getattr(improper.connection_members[2].atom_type, "name", None) + or improper.connection_members[2].name + ) + p4_name = ( + getattr(improper.connection_members[3].atom_type, "name", None) + or improper.connection_members[3].name + ) + improper_connections = [p1_name, p2_name, p3_name, p4_name] + if not tuple(improper_connections[::-1]) in self.improper_types: + self.improper_types.add(tuple(improper_connections)) + + def _identify_topology_information(self, gmso_molecule): + self._identify_particle_information(gmso_molecule) + self._identify_pairs(self.particle_types) + self._identify_bond_types(gmso_molecule) + self._identify_angle_types(gmso_molecule) + self._identify_dihedral_types(gmso_molecule) + self._identify_improper_types(gmso_molecule) + + def _validate_force_field(self): + self.ff_type = None + if isinstance(self.force_field, str): + ff_xml_path, ff_type = find_xml_ff(self.force_field) + self.ff_type = ff_type + self.gmso_molecule = apply_xml_ff(ff_xml_path, self.gmso_molecule) + # Update topology information from typed gmso after applying ff. + self._identify_topology_information(self.gmso_molecule) + elif isinstance(self.force_field, List): + _validate_hoomd_ff(self.force_field, self.topology_information) + self.ff_type = FF_Types.Hoomd + + def assign_mol_name(self, name): + for mol in self.molecules: + mol.name = name + + +class Polymer(Molecule): + def __init__( + self, + lengths, + num_mols, + smiles=None, + file=None, + force_field=None, + bond_indices=None, + bond_length=None, + bond_orientation=None, + **kwargs, + ): + self.lengths = check_return_iterable(lengths) + self.bond_indices = bond_indices + self.bond_length = bond_length + self.bond_orientation = bond_orientation + num_mols = check_return_iterable(num_mols) + if len(num_mols) != len(self.lengths): + raise ValueError("Number of molecules and lengths must be equal.") + super(Polymer, self).__init__( + num_mols=num_mols, + smiles=smiles, + file=file, + force_field=force_field, + **kwargs, + ) + + @property + def monomer(self): + return self._mb_molecule + + def _build(self, length): + chain = mbPolymer() + chain.add_monomer( + self.monomer, + indices=self.bond_indices, + separation=self.bond_length, + orientation=self.bond_orientation, + ) + chain.build(n=length, sequence="A") + return chain + + def _generate(self): + for idx, length in enumerate(self.lengths): + for i in range(self.n_mols[idx]): + mol = self._build(length=length) + self._molecules.append(mol) + + +class CoPolymer(Molecule): + """Builds a polymer consisting of two monomer types. + + Parameters + ---------- + monomer_A : hoomd_polymers.molecules.Polymer; required + Class of the A-type monomer + monomer_B : hoomd_polymers.molecules.Polymer: required + Class of the B-type monomer + length : int; required + The total number of monomers in the molecule + sequence : str; optional; default None + Manually define the sequence of 'A' and 'B' monomers. + Leave as None if generating random sequences. + Example: sequence = "AABAABAAB" + random_sequence : bool; optional; default False + Creates a random 'A' 'B' sequence as a function + of the AB_ratio. + AB_ratio : float; optional; default 0.50 + The relative weight of A to B monomer types. + Used when generating random sequences. + seed : int; optional; default 24 + Set the seed used when generating random sequences + """ + + def __init__( + self, + monomer_A, + monomer_B, + lengths, + num_mols, + force_field=None, + sequence=None, + random_sequence=False, + AB_ratio=0.50, + seed=24, + ): + self.lengths = check_return_iterable(lengths) + self.monomer_A = monomer_A(lengths=[1], num_mols=[1]) + self.monomer_B = monomer_B(lengths=[1], num_mols=[1]) + num_mols = check_return_iterable(num_mols) + if len(num_mols) != len(self.lengths): + raise ValueError("Number of molecules and lengths must be equal.") + self.sequence = sequence + self.random_sequence = random_sequence + self.AB_ratio = AB_ratio + self.seed = seed + self._A_count = 0 + self._B_count = 0 + self.smiles = [self.monomer_A.smiles, self.monomer_B.smiles] + self.file = [self.monomer_A.file, self.monomer_B.file] + random.seed(self.seed) + super(CoPolymer, self).__init__( + num_mols=num_mols, + smiles=self.smiles, + file=self.file, + force_field=force_field, + ) + + @property + def A_ratio(self): + return self._A_count / (self._A_count + self._B_count) + + @property + def B_ratio(self): + return self._B_count / (self._A_count + self._B_count) + + def _build(self, length, sequence): + chain = mbPolymer() + chain.add_monomer( + self.monomer_A.monomer, + indices=self.monomer_A.bond_indices, + orientation=self.monomer_A.bond_orientation, + separation=self.monomer_A.bond_length, + ) + chain.add_monomer( + self.monomer_B.monomer, + indices=self.monomer_B.bond_indices, + orientation=self.monomer_B.bond_orientation, + separation=self.monomer_B.bond_length, + ) + chain.build(n=length, sequence=sequence) + return chain + + def _load(self): + return None + + def _generate(self): + for idx, length in enumerate(self.lengths): + for i in range(self.n_mols[idx]): + if self.random_sequence: + sequence = random.choices( + ["A", "B"], [self.AB_ratio, 1 - self.AB_ratio], k=length + ) + self._A_count += sequence.count("A") + self._B_count += sequence.count("B") + _length = 1 + else: + sequence = self.sequence + _length = length + mol = self._build(length=_length, sequence=sequence) + self._molecules.append(mol) diff --git a/hoomd_polymers/sim/simulation.py b/hoomd_polymers/base/simulation.py similarity index 59% rename from hoomd_polymers/sim/simulation.py rename to hoomd_polymers/base/simulation.py index 7b42f7f4..72e1847b 100644 --- a/hoomd_polymers/sim/simulation.py +++ b/hoomd_polymers/base/simulation.py @@ -1,16 +1,14 @@ -from itertools import combinations_with_replacement as combo -import operator -import os import pickle +import warnings import gsd.hoomd import hoomd import hoomd.md -from mbuild.formats.hoomd_forcefield import create_hoomd_forcefield import numpy as np -import parmed as pmd +import unyt as u -from hoomd_polymers.sim.actions import UpdateWalls, StdOutLogger +from hoomd_polymers.utils import StdOutLogger, UpdateWalls +from hoomd_polymers.utils.exceptions import ReferenceUnitError class Simulation(hoomd.simulation.Simulation): @@ -44,13 +42,12 @@ class Simulation(hoomd.simulation.Simulation): The file name to use for the .txt log file seed : int, default 42 Seed passed to integrator when randomizing velocities. - restart : str, default None - Path to gsd file from which to restart the simulation Methods ------- """ + def __init__( self, initial_state, @@ -62,7 +59,7 @@ def __init__( gsd_write_freq=1e4, gsd_file_name="trajectory.gsd", log_write_freq=1e3, - log_file_name="sim_data.txt" + log_file_name="sim_data.txt", ): super(Simulation, self).__init__(device, seed) self.initial_state = initial_state @@ -70,7 +67,9 @@ def __init__( self.r_cut = r_cut self.gsd_write_freq = int(gsd_write_freq) self.log_write_freq = int(log_write_freq) - self._std_out_freq = int((self.gsd_write_freq + self.log_write_freq)/2) + self._std_out_freq = int( + (self.gsd_write_freq + self.log_write_freq) / 2 + ) self.gsd_file_name = gsd_file_name self.log_file_name = log_file_name self.log_quantities = [ @@ -83,9 +82,7 @@ def __init__( ] self.integrator = None self._dt = dt - self._reference_distance = 1 - self._reference_energy = 1 - self._reference_mass = 1 + self._reference_values = dict() self._integrate_group = hoomd.filter.All() self._wall_forces = dict() self._create_state(self.initial_state) @@ -100,28 +97,79 @@ def forces(self): return self._forcefield @property - def reference_distance(self): - return self._reference_distance + def reference_length(self): + return self._reference_values.get("length", None) - @reference_distance.setter - def reference_distance(self, distance): - self._reference_distance = distance + @property + def reference_mass(self): + return self._reference_values.get("mass", None) @property def reference_energy(self): - return self._reference_energy - - @reference_energy.setter - def reference_energy(self, energy): - self._reference_energy = energy + return self._reference_values.get("energy", None) @property - def reference_mass(self): - return self._reference_mass + def reference_values(self): + return self._reference_values + + @reference_length.setter + def reference_length(self, length, unit=None): + if isinstance(length, u.array.unyt_quantity): + self._reference_values["length"] = length + elif isinstance(unit, str) and ( + isinstance(length, float) or isinstance(length, int) + ): + self._reference_values["length"] = length * getattr(u, unit) + else: + raise ReferenceUnitError( + f"Invalid reference length input.Please provide reference " + f"length (number) and unit (string) or pass length value as an " + f"{str(u.array.unyt_quantity)}." + ) + + @reference_energy.setter + def reference_energy(self, energy, unit=None): + if isinstance(energy, u.array.unyt_quantity): + self._reference_values["energy"] = energy + elif isinstance(unit, str) and ( + isinstance(energy, float) or isinstance(energy, int) + ): + self._reference_values["energy"] = energy * getattr(u, unit) + else: + raise ReferenceUnitError( + f"Invalid reference energy input.Please provide reference " + f"energy (number) and unit (string) or pass energy value as an " + f"{str(u.array.unyt_quantity)}." + ) @reference_mass.setter - def reference_mass(self, mass): - self._reference_mass = mass + def reference_mass(self, mass, unit=None): + if isinstance(mass, u.array.unyt_quantity): + self._reference_values["mass"] = mass + elif isinstance(unit, str) and ( + isinstance(mass, float) or isinstance(mass, int) + ): + self._reference_values["mass"] = mass * getattr(u, unit) + else: + raise ReferenceUnitError( + f"Invalid reference mass input.Please provide reference " + f"mass (number) and " + f"unit (string) or pass mass value as an " + f"{str(u.array.unyt_quantity)}." + ) + + @reference_values.setter + def reference_values(self, ref_value_dict): + ref_keys = ["length", "mass", "energy"] + for k in ref_keys: + if k not in ref_value_dict.keys(): + raise ValueError(f"Missing reference for {k}.") + if not isinstance(ref_value_dict[k], u.array.unyt_quantity): + raise ReferenceUnitError( + f"{k} reference value must be of type " + f"{str(u.array.unyt_quantity)}" + ) + self._reference_values = ref_value_dict @property def box_lengths_reduced(self): @@ -130,7 +178,15 @@ def box_lengths_reduced(self): @property def box_lengths(self): - return self.box_lengths_reduced * self.reference_distance + if self.reference_length: + return self.box_lengths_reduced * self.reference_length + else: + warnings.warn( + "Reference length is not specified. Using HOOMD's unit-less " + "length instead. You can set reference length value and unit " + "with `reference_length()` method. " + ) + return self.box_lengths_reduced @property def volume_reduced(self): @@ -147,15 +203,23 @@ def mass_reduced(self): @property def mass(self): - return self.mass_reduced * self.reference_mass + if self.reference_mass: + return self.mass_reduced * self.reference_mass + else: + warnings.warn( + "Reference mass is not specified. Using HOOMD's unit-less mass " + "instead. You can set reference mass value and unit with " + "`reference_mass()` method. " + ) + return self.mass_reduced @property def density_reduced(self): - return (self.mass_reduced / self.volume_reduced) + return self.mass_reduced / self.volume_reduced @property def density(self): - return (self.mass / self.volume) + return self.mass / self.volume @property def nlist(self): @@ -179,10 +243,19 @@ def dt(self, value): @property def real_timestep(self): - mass = self.reference_mass.to("kg") - dist = self.reference_distance.to("m") - energy = self.reference_energy.to("J") - tau = (mass*(dist**2))/energy + if self._reference_values.get("mass"): + mass = self._reference_values["mass"].to("kg") + else: + mass = 1 * u.kg + if self._reference_values.get("length"): + dist = self.reference_length.to("m") + else: + dist = 1 * u.m + if self._reference_values.get("energy"): + energy = self.reference_energy.to("J") + else: + energy = 1 * u.J + tau = (mass * (dist**2)) / energy timestep = self.dt * (tau**0.5) return timestep @@ -202,9 +275,9 @@ def method(self): return self.operations.integrator.methods[0] else: raise RuntimeError( - "No integrator, or method has been set yet. " - "These will be set once one of the run functions " - "have been called for the first time." + "No integrator, or method has been set yet. " + "These will be set once one of the run functions " + "have been called for the first time." ) def add_force(self, hoomd_force): @@ -225,11 +298,11 @@ def adjust_epsilon(self, scale_by=None, shift_by=None, type_filter=None): for k in lj_forces.params.keys(): if type_filter and k not in type_filter: continue - epsilon = lj_forces.params[k]['epsilon'] + epsilon = lj_forces.params[k]["epsilon"] if scale_by: - lj_forces.params[k]['epsilon'] = epsilon * scale_by + lj_forces.params[k]["epsilon"] = epsilon * scale_by elif shift_by: - lj_forces.params[k]['epsilon'] = epsilon + shift_by + lj_forces.params[k]["epsilon"] = epsilon + shift_by def adjust_sigma(self, scale_by=None, shift_by=None, type_filter=None): """""" @@ -237,11 +310,11 @@ def adjust_sigma(self, scale_by=None, shift_by=None, type_filter=None): for k in lj_forces.params.keys(): if type_filter and k not in type_filter: continue - sigma = lj_forces.params[k]['sigma'] + sigma = lj_forces.params[k]["sigma"] if scale_by: - lj_forces.params[k]['sigma'] = sigma * scale_by + lj_forces.params[k]["sigma"] = sigma * scale_by elif shift_by: - lj_forces.params[k]['sigma'] = sigma + shift_by + lj_forces.params[k]["sigma"] = sigma + shift_by def set_integrator_method(self, integrator_method, method_kwargs): """Creates an initial (or updates the existing) method used by @@ -257,13 +330,13 @@ def set_integrator_method(self, integrator_method, method_kwargs): A diction of parameter:value for the integrator method used """ - if not self.integrator: # Integrator and method not yet created + if not self.integrator: # Integrator and method not yet created self.integrator = hoomd.md.Integrator(dt=self.dt) self.integrator.forces = self._forcefield self.operations.add(self.integrator) new_method = integrator_method(**method_kwargs) self.operations.integrator.methods = [new_method] - else: # Replace the existing integrator method + else: # Replace the existing integrator method self.integrator.methods.remove(self.method) new_method = integrator_method(**method_kwargs) self.integrator.methods.append(new_method) @@ -271,7 +344,7 @@ def set_integrator_method(self, integrator_method, method_kwargs): def add_walls(self, wall_axis, sigma, epsilon, r_cut, r_extrap=0): """""" wall_axis = np.asarray(wall_axis) - wall_origin = wall_axis * self.box_lengths_reduced/2 + wall_origin = wall_axis * self.box_lengths_reduced / 2 wall_normal = -wall_axis wall_origin2 = -wall_origin wall_normal2 = -wall_normal @@ -279,18 +352,20 @@ def add_walls(self, wall_axis, sigma, epsilon, r_cut, r_extrap=0): wall2 = hoomd.wall.Plane(origin=wall_origin2, normal=wall_normal2) lj_walls = hoomd.md.external.wall.LJ(walls=[wall1, wall2]) lj_walls.params[self.state.particle_types] = { - "epsilon": epsilon, - "sigma": sigma, - "r_cut": r_cut, - "r_extrap": r_extrap + "epsilon": epsilon, + "sigma": sigma, + "r_cut": r_cut, + "r_extrap": r_extrap, } self.add_force(lj_walls) self._wall_forces[tuple(wall_axis)] = ( - lj_walls, - {"sigma": sigma, - "epsilon": epsilon, - "r_cut": r_cut, - "r_extrap": r_extrap} + lj_walls, + { + "sigma": sigma, + "epsilon": epsilon, + "r_cut": r_cut, + "r_extrap": r_extrap, + }, ) def remove_walls(self, wall_axis): @@ -299,13 +374,13 @@ def remove_walls(self, wall_axis): self.remove_force(wall_force) def run_update_volume( - self, - n_steps, - period, - kT, - tau_kt, - final_box_lengths, - thermalize_particles=True + self, + n_steps, + period, + kT, + tau_kt, + final_box_lengths, + thermalize_particles=True, ): """Runs an NVT simulation while shrinking or expanding the simulation volume to the given final volume. @@ -326,26 +401,28 @@ def run_update_volume( """ resize_trigger = hoomd.trigger.Periodic(period) box_ramp = hoomd.variant.Ramp( - A=0, B=1, t_start=self.timestep, t_ramp=int(n_steps) + A=0, B=1, t_start=self.timestep, t_ramp=int(n_steps) ) initial_box = self.state.box final_box = hoomd.Box( - Lx=final_box_lengths[0], - Ly=final_box_lengths[1], - Lz=final_box_lengths[2] + Lx=final_box_lengths[0], + Ly=final_box_lengths[1], + Lz=final_box_lengths[2], ) box_resizer = hoomd.update.BoxResize( - box1=initial_box, - box2=final_box, - variant=box_ramp, - trigger=resize_trigger + box1=initial_box, + box2=final_box, + variant=box_ramp, + trigger=resize_trigger, ) self.operations.updaters.append(box_resizer) self.set_integrator_method( - integrator_method=hoomd.md.methods.NVT, - method_kwargs={ - "tau": tau_kt, "filter": self.integrate_group, "kT": kT - }, + integrator_method=hoomd.md.methods.NVT, + method_kwargs={ + "tau": tau_kt, + "filter": self.integrate_group, + "kT": kT, + }, ) if thermalize_particles: self._thermalize_system(kT) @@ -353,86 +430,86 @@ def run_update_volume( if self._wall_forces: wall_update = UpdateWalls(sim=self) wall_updater = hoomd.update.CustomUpdater( - trigger=resize_trigger, action=wall_update + trigger=resize_trigger, action=wall_update ) self.operations.updaters.append(wall_updater) std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( - trigger=hoomd.trigger.Periodic(self._std_out_freq), - action=std_out_logger + trigger=hoomd.trigger.Periodic(self._std_out_freq), + action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) self.run(n_steps + 1) self.operations.updaters.remove(std_out_logger_printer) def run_langevin( - self, - n_steps, - kT, - alpha, - tally_reservoir_energy=False, - default_gamma=1.0, - default_gamma_r=(1.0, 1.0, 1.0), - thermalize_particles=True + self, + n_steps, + kT, + alpha, + tally_reservoir_energy=False, + default_gamma=1.0, + default_gamma_r=(1.0, 1.0, 1.0), + thermalize_particles=True, ): """""" self.set_integrator_method( - integrator_method=hoomd.md.methods.Langevin, - method_kwargs={ - "filter": self.integrate_group, - "kT": kT, - "alpha": alpha, - "tally_reservoir_energy": tally_reservoir_energy, - "default_gamma": default_gamma, - "default_gamma_r": default_gamma_r, - } + integrator_method=hoomd.md.methods.Langevin, + method_kwargs={ + "filter": self.integrate_group, + "kT": kT, + "alpha": alpha, + "tally_reservoir_energy": tally_reservoir_energy, + "default_gamma": default_gamma, + "default_gamma_r": default_gamma_r, + }, ) if thermalize_particles: self._thermalize_system(kT) std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( - trigger=hoomd.trigger.Periodic(self._std_out_freq), - action=std_out_logger + trigger=hoomd.trigger.Periodic(self._std_out_freq), + action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) self.run(n_steps) self.operations.updaters.remove(std_out_logger_printer) def run_NPT( - self, - n_steps, - kT, - pressure, - tau_kt, - tau_pressure, - couple="xyz", - box_dof=[True, True, True, False, False, False], - rescale_all=False, - gamma=0.0, - thermalize_particles=True + self, + n_steps, + kT, + pressure, + tau_kt, + tau_pressure, + couple="xyz", + box_dof=[True, True, True, False, False, False], + rescale_all=False, + gamma=0.0, + thermalize_particles=True, ): """""" self.set_integrator_method( - integrator_method=hoomd.md.methods.NPT, - method_kwargs={ - "kT": kT, - "S": pressure, - "tau": tau_kt, - "tauS": tau_pressure, - "couple": couple, - "box_dof": box_dof, - "rescale_all": rescale_all, - "gamma": gamma, - "filter": self.integrate_group, - "kT": kT - } + integrator_method=hoomd.md.methods.NPT, + method_kwargs={ + "kT": kT, + "S": pressure, + "tau": tau_kt, + "tauS": tau_pressure, + "couple": couple, + "box_dof": box_dof, + "rescale_all": rescale_all, + "gamma": gamma, + "filter": self.integrate_group, + "kT": kT, + }, ) if thermalize_particles: self._thermalize_system(kT) std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( - trigger=hoomd.trigger.Periodic(self._std_out_freq), - action=std_out_logger + trigger=hoomd.trigger.Periodic(self._std_out_freq), + action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) self.run(n_steps) @@ -441,17 +518,19 @@ def run_NPT( def run_NVT(self, n_steps, kT, tau_kt, thermalize_particles=True): """""" self.set_integrator_method( - integrator_method=hoomd.md.methods.NVT, - method_kwargs={ - "tau": tau_kt, "filter": self.integrate_group, "kT": kT - } + integrator_method=hoomd.md.methods.NVT, + method_kwargs={ + "tau": tau_kt, + "filter": self.integrate_group, + "kT": kT, + }, ) if thermalize_particles: self._thermalize_system(kT) std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( - trigger=hoomd.trigger.Periodic(self._std_out_freq), - action=std_out_logger + trigger=hoomd.trigger.Periodic(self._std_out_freq), + action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) self.run(n_steps) @@ -460,25 +539,28 @@ def run_NVT(self, n_steps, kT, tau_kt, thermalize_particles=True): def run_NVE(self, n_steps): """""" self.set_integrator_method( - integrator_method=hoomd.md.methods.NVE, - method_kwargs={"filter": self.integrate_group} + integrator_method=hoomd.md.methods.NVE, + method_kwargs={"filter": self.integrate_group}, ) std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( - trigger=hoomd.trigger.Periodic(self._std_out_freq), - action=std_out_logger + trigger=hoomd.trigger.Periodic(self._std_out_freq), + action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) self.run(n_steps) self.operations.updaters.remove(std_out_logger_printer) def run_displacement_cap(self, n_steps, maximum_displacement=1e-3): - """ NVE based integrator that Puts a cap on the maximum displacement per time step. - - DisplacementCapped method is mostly useful for initially relaxing a system with overlapping particles. - Putting a cap on the max particle displacement prevents Hoomd Particle Out of Box execption. - Once the system is relaxed, other run methods (NVE, NVT, etc) can be used. - + """NVE based integrator that Puts a cap on the maximum displacement + per time step. + + DisplacementCapped method is mostly useful for initially relaxing a + system with overlapping particles. Putting a cap on the max particle + displacement prevents Hoomd Particle Out of Box execption. + Once the system is relaxed, other run methods (NVE, NVT, etc) can be + used. + Parameters: ----------- n_steps : int, required @@ -487,14 +569,16 @@ def run_displacement_cap(self, n_steps, maximum_displacement=1e-3): """ self.set_integrator_method( - integrator_method=hoomd.md.methods.DisplacementCapped, - method_kwargs={"filter": self.integrate_group, - "maximum_displacement": maximum_displacement} + integrator_method=hoomd.md.methods.DisplacementCapped, + method_kwargs={ + "filter": self.integrate_group, + "maximum_displacement": maximum_displacement, + }, ) std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( - trigger=hoomd.trigger.Periodic(self._std_out_freq), - action=std_out_logger + trigger=hoomd.trigger.Periodic(self._std_out_freq), + action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) self.run(n_steps) @@ -502,10 +586,7 @@ def run_displacement_cap(self, n_steps, maximum_displacement=1e-3): def temperature_ramp(self, n_steps, kT_start, kT_final): return hoomd.variant.Ramp( - A=kT_start, - B=kT_final, - t_start=self.timestep, - t_ramp=int(n_steps) + A=kT_start, B=kT_final, t_start=self.timestep, t_ramp=int(n_steps) ) def pickle_forcefield(self, file_path="forcefield.pickle"): @@ -518,26 +599,30 @@ def save_restart_gsd(self, file_path="restart.gsd"): def _thermalize_system(self, kT): if isinstance(kT, hoomd.variant.Ramp): self.state.thermalize_particle_momenta( - filter=self.integrate_group, kT=kT.range[0] + filter=self.integrate_group, kT=kT.range[0] ) else: self.state.thermalize_particle_momenta( - filter=self.integrate_group, kT=kT + filter=self.integrate_group, kT=kT ) def _lj_force(self): if not self.integrator: lj_force = [ - f for f in self._forcefield if - isinstance(f, hoomd.md.pair.pair.LJ)][0] + f + for f in self._forcefield + if isinstance(f, hoomd.md.pair.pair.LJ) + ][0] else: lj_force = [ - f for f in self.integrator.forces if - isinstance(f, hoomd.md.pair.pair.LJ)][0] + f + for f in self.integrator.forces + if isinstance(f, hoomd.md.pair.pair.LJ) + ][0] return lj_force def _create_state(self, initial_state): - if isinstance(initial_state, str): # Load from a GSD file + if isinstance(initial_state, str): # Load from a GSD file print("Initializing simulation state from a GSD file.") self.create_state_from_gsd(initial_state) elif isinstance(initial_state, hoomd.snapshot.Snapshot): @@ -550,16 +635,16 @@ def _create_state(self, initial_state): def _add_hoomd_writers(self): """Creates gsd and log writers""" gsd_writer = hoomd.write.GSD( - filename=self.gsd_file_name, - trigger=hoomd.trigger.Periodic(int(self.gsd_write_freq)), - mode="wb", - dynamic=["momentum"] + filename=self.gsd_file_name, + trigger=hoomd.trigger.Periodic(int(self.gsd_write_freq)), + mode="wb", + dynamic=["momentum"], ) logger = hoomd.logging.Logger(categories=["scalar", "string"]) logger.add(self, quantities=["timestep", "tps"]) thermo_props = hoomd.md.compute.ThermodynamicQuantities( - filter=self.integrate_group + filter=self.integrate_group ) self.operations.computes.append(thermo_props) logger.add(thermo_props, quantities=self.log_quantities) diff --git a/hoomd_polymers/base/system.py b/hoomd_polymers/base/system.py new file mode 100644 index 00000000..0a8fe8ee --- /dev/null +++ b/hoomd_polymers/base/system.py @@ -0,0 +1,526 @@ +import warnings +from abc import ABC, abstractmethod +from typing import List + +import gsd +import mbuild as mb +import numpy as np +import unyt as u +from gmso.external import ( + from_mbuild, + from_parmed, + to_gsd_snapshot, + to_hoomd_forcefield, + to_parmed, +) +from gmso.parameterization import apply + +from hoomd_polymers.base.molecule import Molecule +from hoomd_polymers.utils import FF_Types, check_return_iterable, xml_to_gmso_ff +from hoomd_polymers.utils.exceptions import ( + ForceFieldError, + MoleculeLoadError, + ReferenceUnitError, +) + + +class System(ABC): + """Base class from which other systems inherit. + + Parameters + ---------- + molecule : hoomd_polymers.molecule; required + n_mols : int; required + The number of times to replicate molecule in the system + density : float; optional; default None + The desired density of the system (g/cm^3). Used to set the + target_box attribute. Can be useful when initializing + systems at low denisty and running a shrink simulaton + to acheive a target density. + """ + + def __init__( + self, + molecules, + density: float, + r_cut: float, + force_field=None, + auto_scale=False, + remove_hydrogens=False, + remove_charges=False, + scale_charges=False, + base_units=dict(), + ): + self._molecules = check_return_iterable(molecules) + self._force_field = None + if force_field: + self._force_field = check_return_iterable(force_field) + self.density = density + self.r_cut = r_cut + self.auto_scale = auto_scale + self.remove_hydrogens = remove_hydrogens + self.remove_charges = remove_charges + self.scale_charges = scale_charges + self.target_box = None + self.all_molecules = [] + self._hoomd_snapshot = None + self._hoomd_forcefield = [] + self._reference_values = base_units + self._gmso_forcefields_dict = dict() + self.gmso_system = None + + # Collecting all molecules + self.n_mol_types = 0 + for mol_item in self._molecules: + if isinstance(mol_item, Molecule): + mol_item.assign_mol_name(str(self.n_mol_types)) + self.all_molecules.extend(mol_item.molecules) + # if ff is provided in Molecule class + if mol_item.force_field: + if mol_item.ff_type == FF_Types.Hoomd: + self._hoomd_forcefield.extend(mol_item.force_field) + else: + self._gmso_forcefields_dict[ + str(self.n_mol_types) + ] = xml_to_gmso_ff(mol_item.force_field) + self.n_mol_types += 1 + elif isinstance(mol_item, mb.Compound): + mol_item.name = str(self.n_mol_types) + self.all_molecules.append(mol_item) + elif isinstance(mol_item, List): + for sub_mol in mol_item: + if isinstance(sub_mol, mb.Compound): + sub_mol.name = str(self.n_mol_types) + self.all_molecules.append(sub_mol) + else: + raise MoleculeLoadError( + msg=f"Unsupported compound type {type(sub_mol)}. " + f"Supported compound types are: " + f"{str(mb.Compound)}" + ) + self.n_mol_types += 1 + + # Collecting all force-fields only if xml force-field is provided + if self._force_field: + for i in range(self.n_mol_types): + if not self._gmso_forcefields_dict.get(str(i)): + if i < len(self._force_field): + # if there is a ff for each molecule type + ff_index = i + else: + # if there is only one ff for all molecule types + ff_index = 0 + if getattr(self._force_field[ff_index], "gmso_ff"): + self._gmso_forcefields_dict[str(i)] = self._force_field[ + ff_index + ].gmso_ff + else: + raise ForceFieldError( + msg=f"GMSO Force field in " + f"{self._force_field[ff_index]} is not " + f"provided." + ) + self.system = self._build_system() + self.gmso_system = self._convert_to_gmso() + if self._force_field: + self._apply_forcefield() + if self.remove_hydrogens: + self._remove_hydrogens() + + @abstractmethod + def _build_system(self): + pass + + @property + def n_molecules(self): + return len(self.all_molecules) + + @property + def n_particles(self): + if self.gmso_system: + return self.gmso_system.n_sites + return sum([mol.n_particles for mol in self.all_molecules]) + + @property + def mass(self): + if self.gmso_system: + return sum( + float(site.mass.to("amu").value) + for site in self.gmso_system.sites + ) + return sum(mol.mass for mol in self.all_molecules) + + @property + def box(self): + return self.system.box + + @property + def reference_length(self): + return self._reference_values.get("length", None) + + @property + def reference_mass(self): + return self._reference_values.get("mass", None) + + @property + def reference_energy(self): + return self._reference_values.get("energy", None) + + @property + def reference_values(self): + return self._reference_values + + @reference_length.setter + def reference_length(self, length, unit=None): + if isinstance(length, u.array.unyt_quantity): + self._reference_values["length"] = length + elif isinstance(unit, str) and ( + isinstance(length, float) or isinstance(length, int) + ): + self._reference_values["length"] = length * getattr(u, unit) + else: + raise ReferenceUnitError( + f"Invalid reference length input.Please provide reference " + f"length (number) and unit (string) or pass length value as " + f"an {str(u.array.unyt_quantity)}." + ) + + @reference_energy.setter + def reference_energy(self, energy, unit=None): + if isinstance(energy, u.array.unyt_quantity): + self._reference_values["energy"] = energy + elif isinstance(unit, str) and ( + isinstance(energy, float) or isinstance(energy, int) + ): + self._reference_values["energy"] = energy * getattr(u, unit) + else: + raise ReferenceUnitError( + f"Invalid reference energy input.Please provide reference " + f"energy (number) and unit (string) or pass energy value as " + f"an {str(u.array.unyt_quantity)}." + ) + + @reference_mass.setter + def reference_mass(self, mass, unit=None): + if isinstance(mass, u.array.unyt_quantity): + self._reference_values["mass"] = mass + elif isinstance(unit, str) and ( + isinstance(mass, float) or isinstance(mass, int) + ): + self._reference_values["mass"] = mass * getattr(u, unit) + else: + raise ReferenceUnitError( + f"Invalid reference mass input.Please provide reference " + f"mass (number) and unit (string) or pass mass value as an " + f"{str(u.array.unyt_quantity)}." + ) + + @reference_values.setter + def reference_values(self, ref_value_dict): + ref_keys = ["length", "mass", "energy"] + for k in ref_keys: + if k not in ref_value_dict.keys(): + raise ValueError(f"Missing reference for {k}.") + if not isinstance(ref_value_dict[k], u.array.unyt_quantity): + raise ReferenceUnitError( + f"{k} reference value must be of type " + f"{str(u.array.unyt_quantity)}" + ) + self._reference_values = ref_value_dict + + @property + def hoomd_snapshot(self): + self._hoomd_snapshot = self._create_hoomd_snapshot() + return self._hoomd_snapshot + + @property + def hoomd_forcefield(self): + self._hoomd_forcefield = self._create_hoomd_forcefield() + return self._hoomd_forcefield + + def _remove_hydrogens(self): + """Call this method to remove hydrogen atoms from the system. + The masses and charges of the hydrogens are absorbed into + the heavy atoms they were bonded to. + """ + parmed_struc = to_parmed(self.gmso_system) + # Try by element first: + hydrogens = [a for a in parmed_struc.atoms if a.element == 1] + if len(hydrogens) == 0: # Try by mass + hydrogens = [a for a in parmed_struc.atoms if a.mass == 1.008] + if len(hydrogens) == 0: + warnings.warn( + "Hydrogen atoms could not be found by element or mass" + ) + for h in hydrogens: + h.atomic_number = 1 + bonded_atom = h.bond_partners[0] + bonded_atom.mass += h.mass + bonded_atom.charge += h.charge + parmed_struc.strip([a.atomic_number == 1 for a in parmed_struc.atoms]) + if len(hydrogens) > 0: + self.gmso_system = from_parmed(parmed_struc) + + def to_gsd(self, file_name): + with gsd.hoomd.open(file_name, "wb") as traj: + traj.append(self.hoomd_snapshot) + + def _convert_to_gmso(self): + topology = from_mbuild(self.system) + topology.identify_connections() + return topology + + def _create_hoomd_forcefield(self): + force_list = [] + ff, refs = to_hoomd_forcefield( + top=self.gmso_system, + r_cut=self.r_cut, + auto_scale=self.auto_scale, + base_units=self._reference_values + if self._reference_values + else None, + ) + for force in ff: + force_list.extend(ff[force]) + return force_list + + def _create_hoomd_snapshot(self): + snap, refs = to_gsd_snapshot( + top=self.gmso_system, + auto_scale=self.auto_scale, + base_units=self._reference_values + if self._reference_values + else None, + ) + return snap + + def _apply_forcefield(self): + self.gmso_system = apply( + self.gmso_system, + self._gmso_forcefields_dict, + identify_connections=True, + use_molecule_info=True, + identify_connected_components=False, + ) + if self.remove_charges: + for site in self.gmso_system.sites: + site.charge = 0 + if self.scale_charges and not self.remove_charges: + pass + # TODO: Scale charges from self.gmso_system + epsilons = [ + s.atom_type.parameters["epsilon"] for s in self.gmso_system.sites + ] + sigmas = [ + s.atom_type.parameters["sigma"] for s in self.gmso_system.sites + ] + masses = [s.mass for s in self.gmso_system.sites] + + energy_scale = np.max(epsilons) if self.auto_scale else 1.0 + length_scale = np.max(sigmas) if self.auto_scale else 1.0 + mass_scale = np.max(masses) if self.auto_scale else 1.0 + + self._reference_values["energy"] = energy_scale * epsilons[0].unit_array + self._reference_values["length"] = length_scale * sigmas[0].unit_array + if self.auto_scale: + self._reference_values["mass"] = mass_scale * masses[ + 0 + ].unit_array.to("amu") + else: + self._reference_values["mass"] = mass_scale * u.g / u.mol + + def set_target_box( + self, x_constraint=None, y_constraint=None, z_constraint=None + ): + """Set the target volume of the system during + the initial shrink step. + If no constraints are set, the target box is cubic. + Setting constraints will hold those box vectors + constant and adjust others to match the target density. + + Parameters + ----------- + x_constraint : float, optional, defualt=None + Fixes the box length (nm) along the x axis + y_constraint : float, optional, default=None + Fixes the box length (nm) along the y axis + z_constraint : float, optional, default=None + Fixes the box length (nm) along the z axis + + """ + if not any([x_constraint, y_constraint, z_constraint]): + Lx = Ly = Lz = self._calculate_L() + else: + constraints = np.array([x_constraint, y_constraint, z_constraint]) + fixed_L = constraints[np.where(constraints is not None)] + # Conv from nm to cm for _calculate_L + fixed_L *= 1e-7 + L = self._calculate_L(fixed_L=fixed_L) + constraints[np.where(constraints is None)] = L + Lx, Ly, Lz = constraints + + self.target_box = np.array([Lx, Ly, Lz]) + + def visualize(self): + if self.system: + self.system.visualize().show() + else: + raise ValueError( + "The initial configuraiton has not been created yet." + ) + + def _calculate_L(self, fixed_L=None): + """Calculates the required box length(s) given the + mass of a sytem and the target density. + + Box edge length constraints can be set by set_target_box(). + If constraints are set, this will solve for the required + lengths of the remaining non-constrained edges to match + the target density. + + Parameters + ---------- + fixed_L : np.array, optional, defualt=None + Array of fixed box lengths to be accounted for + when solving for L + + """ + # Convert from amu to grams + M = self.mass * 1.66054e-24 + vol = M / self.density # cm^3 + if fixed_L is None: + L = vol ** (1 / 3) + else: + L = vol / np.prod(fixed_L) + if len(fixed_L) == 1: # L is cm^2 + L = L ** (1 / 2) + # Convert from cm back to nm + L *= 1e7 + return L + + +class Pack(System): + """Uses PACKMOL via mbuild.packing.fill_box. + The box used for packing is expanded to allow PACKMOL + to more easily place all the molecules. + + Parameters + ---------- + packing_expand_factor : int; optional, default 5 + + """ + + def __init__( + self, + molecules, + density: float, + r_cut: float, + force_field=None, + auto_scale=False, + remove_hydrogens=False, + remove_charges=False, + scale_charges=False, + base_units=dict(), + packing_expand_factor=5, + edge=0.2, + ): + self.packing_expand_factor = packing_expand_factor + self.edge = edge + super(Pack, self).__init__( + molecules=molecules, + density=density, + force_field=force_field, + r_cut=r_cut, + auto_scale=auto_scale, + remove_hydrogens=remove_hydrogens, + remove_charges=remove_charges, + scale_charges=scale_charges, + base_units=base_units, + ) + + def _build_system(self): + self.set_target_box() + system = mb.packing.fill_box( + compound=self.all_molecules, + n_compounds=[1 for i in self.all_molecules], + box=list(self.target_box * self.packing_expand_factor), + overlap=0.2, + edge=self.edge, + ) + return system + + +class Lattice(System): + """Places the molecules in a lattice configuration. + Assumes two molecules per unit cell. + + Parameters + ---------- + x : float; required + The distance (nm) between lattice points in the x direction. + y : float; required + The distance (nm) between lattice points in the y direction. + n : int; required + The number of times to repeat the unit cell in x and y + lattice_vector : array-like + The vector between points in the unit cell + """ + + def __init__( + self, + molecules, + density: float, + r_cut: float, + x: float, + y: float, + n: int, + basis_vector=[0.5, 0.5, 0], + force_field=None, + auto_scale=False, + remove_hydrogens=False, + remove_charges=False, + scale_charges=False, + base_units=dict(), + ): + self.x = x + self.y = y + self.n = n + self.basis_vector = basis_vector + super(Lattice, self).__init__( + molecules=molecules, + density=density, + force_field=force_field, + r_cut=r_cut, + auto_scale=auto_scale, + remove_hydrogens=remove_hydrogens, + remove_charges=remove_charges, + scale_charges=scale_charges, + base_units=base_units, + ) + + def _build_system(self): + next_idx = 0 + system = mb.Compound() + for i in range(self.n): + layer = mb.Compound() + for j in range(self.n): + try: + comp1 = self.all_molecules[next_idx] + comp2 = self.all_molecules[next_idx + 1] + comp2.translate(self.basis_vector) + # TODO: what if comp1 and comp2 have different names? + unit_cell = mb.Compound( + subcompounds=[comp1, comp2], name=comp1.name + ) + unit_cell.translate((0, self.y * j, 0)) + layer.add(unit_cell) + next_idx += 2 + except IndexError: + pass + layer.translate((self.x * i, 0, 0)) + system.add(layer) + bounding_box = system.get_boundingbox() + x_len = bounding_box.lengths[0] + y_len = bounding_box.lengths[1] + self.set_target_box(x_constraint=x_len, y_constraint=y_len) + return system diff --git a/hoomd_polymers/forcefields.py b/hoomd_polymers/forcefields.py deleted file mode 100644 index 2453a959..00000000 --- a/hoomd_polymers/forcefields.py +++ /dev/null @@ -1,39 +0,0 @@ -import foyer -from hoomd_polymers.library import FF_DIR - - -class GAFF(foyer.Forcefield): - def __init__(self, forcefield_files=f"{FF_DIR}/gaff.xml"): - super(GAFF, self).__init__(forcefield_files=forcefield_files) - self.description = ( - "The General Amber Forcefield written in foyer XML format. " - "The XML file was obtained from the antefoyer package: " - "https://github.com/rsdefever/antefoyer/tree/master/antefoyer" - ) - - -class OPLS_AA(foyer.Forcefield): - def __init__(self, name="oplsaa"): - super(OPLS_AA, self).__init__(name=name) - self.description = ( - "opls-aa forcefield found in the Foyer package." - ) - - -class OPLS_AA_PPS(foyer.Forcefield): - def __init__(self, forcefield_files=f"{FF_DIR}/pps_opls.xml"): - super(OPLS_AA_PPS, self).__init__(forcefield_files=forcefield_files) - self.description = ( - "Based on hoomd_polymers.forcefields.OPLS_AA. " - "Trimmed down to include only PPS parameters. " - "One missing parameter was added manually: " - " " - "The equilibrium angle was determined from " - "experimental PPS papers. The spring constant taken " - "from the equivalent angle in GAFF." - ) - - -class FF_from_file(foyer.Forcefield): - def __init__(self, xml_file): - super(FF_from_file, self).__init__(forcefield_files=xml_file) diff --git a/hoomd_polymers/library/__init__.py b/hoomd_polymers/library/__init__.py index fef490a9..a98cc3cb 100644 --- a/hoomd_polymers/library/__init__.py +++ b/hoomd_polymers/library/__init__.py @@ -1,2 +1,19 @@ -from hoomd_polymers.library.forcefields import FF_DIR -from hoomd_polymers.library.monomers import MON_DIR +from .forcefields import ( + GAFF, + OPLS_AA, + OPLS_AA_BENZENE, + OPLS_AA_DIMETHYLETHER, + OPLS_AA_PPS, + BeadSpring, + FF_from_file, +) +from .polymers import ( + PEEK, + PEKK, + PPS, + LJChain, + PEKK_meta, + PEKK_para, + PolyEthylene, +) +from .simulations.tensile import Tensile diff --git a/hoomd_polymers/library/forcefields.py b/hoomd_polymers/library/forcefields.py new file mode 100644 index 00000000..2e811bd3 --- /dev/null +++ b/hoomd_polymers/library/forcefields.py @@ -0,0 +1,126 @@ +import itertools + +import forcefield_utilities as ffutils +import foyer +import hoomd + +from hoomd_polymers.assets import FF_DIR + + +class GAFF(foyer.Forcefield): + def __init__(self, forcefield_files=f"{FF_DIR}/gaff.xml"): + super(GAFF, self).__init__(forcefield_files=forcefield_files) + self.description = ( + "The General Amber Forcefield written in foyer XML format. " + "The XML file was obtained from the antefoyer package: " + "https://github.com/rsdefever/antefoyer/tree/master/antefoyer" + ) + self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() + + +class OPLS_AA(foyer.Forcefield): + def __init__(self, name="oplsaa"): + super(OPLS_AA, self).__init__(name=name) + self.description = "opls-aa forcefield found in the Foyer package." + self.gmso_ff = ffutils.FoyerFFs().load(name).to_gmso_ff() + + +class OPLS_AA_PPS(foyer.Forcefield): + def __init__(self, forcefield_files=f"{FF_DIR}/pps_opls.xml"): + super(OPLS_AA_PPS, self).__init__(forcefield_files=forcefield_files) + self.description = ( + "Based on hoomd_polymers.forcefields.OPLS_AA. " + "Trimmed down to include only PPS parameters. " + "One missing parameter was added manually: " + " " + "The equilibrium angle was determined from " + "experimental PPS papers. The spring constant taken " + "from the equivalent angle in GAFF." + ) + self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() + + +class OPLS_AA_BENZENE(foyer.Forcefield): + def __init__(self, forcefield_files=f"{FF_DIR}/benzene_opls.xml"): + super(OPLS_AA_BENZENE, self).__init__(forcefield_files=forcefield_files) + self.description = ( + "Based on hoomd_polymers.forcefields.OPLS_AA. " + "Trimmed down to include only benzene parameters." + ) + self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() + + +class OPLS_AA_DIMETHYLETHER(foyer.Forcefield): + def __init__(self, forcefield_files=f"{FF_DIR}/dimethylether_opls.xml"): + super(OPLS_AA_DIMETHYLETHER, self).__init__( + forcefield_files=forcefield_files + ) + self.description = ( + "Based on hoomd_polymers.forcefields.OPLS_AA. " + "Trimmed down to include only dimethyl ether parameters." + ) + self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() + + +class FF_from_file(foyer.Forcefield): + def __init__(self, xml_file): + super(FF_from_file, self).__init__(forcefield_files=xml_file) + self.gmso_ff = ffutils.FoyerFFs().load(xml_file).to_gmso_ff() + + +class BeadSpring: + def __init__( + self, + r_cut, + beads, + bonds=None, + angles=None, + dihedrals=None, + exclusions=["bond", "1-3"], + ): + self.beads = beads + self.bonds = bonds + self.angles = angles + self.dihedrals = dihedrals + self.r_cut = r_cut + self.exclusions = exclusions + self.hoomd_forcefield = self._create_forcefield() + + def _create_forcefield(self): + forces = [] + # Create pair force: + nlist = hoomd.md.nlist.Cell(buffer=0.40, exclusions=self.exclusions) + lj = hoomd.md.pair.LJ(nlist=nlist) + bead_types = [key for key in self.beads.keys()] + all_pairs = list(itertools.combinations_with_replacement(bead_types, 2)) + for pair in all_pairs: + epsilon0 = self.beads[pair[0]]["epsilon"] + epsilon1 = self.beads[pair[1]]["epsilon"] + pair_epsilon = (epsilon0 + epsilon1) / 2 + + sigma0 = self.beads[pair[0]]["sigma"] + sigma1 = self.beads[pair[1]]["sigma"] + pair_sigma = (sigma0 + sigma1) / 2 + + lj.params[pair] = dict(epsilon=pair_epsilon, sigma=pair_sigma) + lj.r_cut[pair] = self.r_cut + forces.append(lj) + # Create bond-stretching force: + if self.bonds: + harmonic_bond = hoomd.md.bond.Harmonic() + for bond_type in self.bonds: + harmonic_bond.params[bond_type] = self.bonds[bond_type] + forces.append(harmonic_bond) + # Create bond-bending force: + if self.angles: + harmonic_angle = hoomd.md.angle.Harmonic() + for angle_type in self.angles: + harmonic_angle.params[angle_type] = self.angles[angle_type] + forces.append(harmonic_angle) + # Create torsion force: + if self.dihedrals: + periodic_dihedral = hoomd.md.dihedral.Periodic() + for dih_type in self.dihedrals: + periodic_dihedral.params[dih_type] = self.dihedrals[dih_type] + forces.append(periodic_dihedral) + return forces diff --git a/hoomd_polymers/library/polymers.py b/hoomd_polymers/library/polymers.py new file mode 100644 index 00000000..3e676803 --- /dev/null +++ b/hoomd_polymers/library/polymers.py @@ -0,0 +1,219 @@ +import os + +import mbuild as mb +from mbuild.coordinate_transform import z_axis_transform + +from hoomd_polymers import CoPolymer, Polymer +from hoomd_polymers.assets import MON_DIR + + +class PolyEthylene(Polymer): + """Creates a Poly(ethylene) chain. + + Parameters + ---------- + length : int; required + The number of monomer repeat units in the chain + """ + + def __init__(self, lengths, num_mols, **kwargs): + smiles = "CC" + bond_indices = [2, 6] + bond_length = 0.145 + bond_orientation = [None, None] + super(PolyEthylene, self).__init__( + lengths=lengths, + num_mols=num_mols, + smiles=smiles, + bond_indices=bond_indices, + bond_length=bond_length, + bond_orientation=bond_orientation, + **kwargs, + ) + + +class PPS(Polymer): + """Creates a Poly(phenylene-sulfide) (PPS) chain. + + Parameters + ---------- + length : int; required + The number of monomer repeat units in the chain + """ + + def __init__(self, lengths, num_mols, **kwargs): + smiles = "c1ccc(S)cc1" + file = None + bond_indices = [7, 10] + bond_length = 0.176 + bond_orientation = [[0, 0, 1], [0, 0, -1]] + super(PPS, self).__init__( + lengths=lengths, + num_mols=num_mols, + smiles=smiles, + file=file, + bond_indices=bond_indices, + bond_length=bond_length, + bond_orientation=bond_orientation, + **kwargs, + ) + + def _load(self): + monomer = mb.load(self.smiles, smiles=True) + # Need to align monomer along zx plane due to orientation of S-H bond + z_axis_transform( + monomer, point_on_z_axis=monomer[7], point_on_zx_plane=monomer[4] + ) + return monomer + + +class PEEK(Polymer): + def __init__(self, lengths, num_mols, **kwargs): + super(PEEK, self).__init__( + lengths=lengths, + num_mols=num_mols, + ) + + +class PEKK(CoPolymer): + def __init__( + self, + lengths, + num_mols, + force_field=None, + sequence=None, + random_sequence=False, + TI_ratio=0.50, + seed=24, + **kwargs, + ): + super(PEKK, self).__init__( + monomer_A=PEKK_meta, + monomer_B=PEKK_para, + lengths=lengths, + num_mols=num_mols, + force_field=force_field, + sequence=sequence, + random_sequence=random_sequence, + AB_ratio=TI_ratio, + seed=seed, + **kwargs, + ) + + +class PEKK_para(Polymer): + """Creates a Poly(ether-ketone-ketone) (PEKK) chain. + The bonding positions of consecutive ketone groups + takes place on the para site of the phenyl ring. + + Parameters + ---------- + length : int; required + The number of monomer repeat units in the chain + """ + + def __init__(self, lengths, num_mols): + smiles = "c1ccc(Oc2ccc(C(=O)c3ccc(C(=O))cc3)cc2)cc1" + file = os.path.join(MON_DIR, "pekk_para.mol2") + bond_indices = [35, 36] + bond_length = 0.148 + bond_orientation = [[0, 0, -1], [0, 0, 1]] + super(PEKK_para, self).__init__( + lengths=lengths, + num_mols=num_mols, + smiles=smiles, + file=file, + bond_indices=bond_indices, + bond_length=bond_length, + bond_orientation=bond_orientation, + ) + + +class PEKK_meta(Polymer): + """Creates a Poly(ether-ketone-ketone) (PEKK) chain. + The bonding positions of consecutive ketone groups + takes place on the meta site of the phenyl ring. + + Parameters + ---------- + length : int; required + The number of monomer repeat units in the chain + """ + + def __init__(self, lengths, num_mols): + smiles = "c1cc(Oc2ccc(C(=O)c3cc(C(=O))ccc3)cc2)ccc1" + file = os.path.join(MON_DIR, "pekk_meta.mol2") + bond_indices = [35, 36] + bond_length = 0.148 + bond_orientation = [[0, 0, -1], [0, 0, 1]] + super(PEKK_meta, self).__init__( + lengths=lengths, + num_mols=num_mols, + smiles=smiles, + file=file, + bond_indices=bond_indices, + bond_length=bond_length, + bond_orientation=bond_orientation, + ) + + +class LJChain(Polymer): + """Creates a coarse-grained bead-spring polymer chain. + + Parameters + ---------- + length : int; required + The number of times to repeat bead_sequence in a single chain. + bead_sequence : list; optional; default ["A"] + The sequence of bead types in the chain. + bond_length : dict; optional; default {"A-A": 1.0} + The bond length between connected beads (units: nm) + bead_mass : dict; optional; default {"A": 1.0} + The mass of the bead types + """ + + def __init__( + self, + lengths, + num_mols, + bead_sequence=["A"], + bead_mass={"A": 1.0}, + bond_lengths={"A-A": 1.0}, + ): + self.bead_sequence = bead_sequence + self.bead_mass = bead_mass + self.bond_lengths = bond_lengths + super(LJChain, self).__init__(lengths=lengths, num_mols=num_mols) + + def _build(self, length): + chain = mb.Compound() + last_bead = None + for i in range(length): + for idx, bead_type in enumerate(self.bead_sequence): + mass = self.bead_mass.get(bead_type, None) + if not mass: + raise ValueError( + f"The bead mass for {bead_type} was not given " + "in the bead_mass dict." + ) + next_bead = mb.Compound(mass=mass, name=bead_type, charge=0) + chain.add(next_bead) + if last_bead: + bead_pair = "-".join([last_bead.name, next_bead.name]) + bond_length = self.bond_lengths.get(bead_pair, None) + if not bond_length: + bead_pair_rev = "-".join( + [next_bead.name, last_bead.name] + ) + bond_length = self.bond_lengths.get(bead_pair_rev, None) + if not bond_length: + raise ValueError( + "The bond length for pair " + f"{bead_pair} or {bead_pair_rev} " + "is not found in the bond_lengths dict." + ) + new_pos = last_bead.xyz[0] + (0, 0, bond_length) + next_bead.translate_to(new_pos) + chain.add_bond([next_bead, last_bead]) + last_bead = next_bead + return chain diff --git a/hoomd_polymers/sim/tensile.py b/hoomd_polymers/library/simulations/tensile.py similarity index 53% rename from hoomd_polymers/sim/tensile.py rename to hoomd_polymers/library/simulations/tensile.py index 3b497f0b..fd63a648 100644 --- a/hoomd_polymers/sim/tensile.py +++ b/hoomd_polymers/library/simulations/tensile.py @@ -1,44 +1,45 @@ import hoomd import numpy as np -from hoomd_polymers.sim.simulation import Simulation -from hoomd_polymers.sim.actions import PullParticles +from hoomd_polymers.base.simulation import Simulation +from hoomd_polymers.utils.actions import PullParticles class Tensile(Simulation): def __init__( - self, - initial_state, - forcefield, - tensile_axis, - fix_ratio=0.20, - r_cut=2.5, - dt=0.0001, - device=hoomd.device.auto_select(), - seed=42, - restart=None, - gsd_write_freq=1e4, - gsd_file_name="trajectory.gsd", - log_write_freq=1e3, - log_file_name="sim_data.txt" + self, + initial_state, + forcefield, + tensile_axis, + fix_ratio=0.20, + r_cut=2.5, + dt=0.0001, + device=hoomd.device.auto_select(), + seed=42, + restart=None, + gsd_write_freq=1e4, + gsd_file_name="trajectory.gsd", + log_write_freq=1e3, + log_file_name="sim_data.txt", ): super(Tensile, self).__init__( - initial_state=initial_state, - forcefield=forcefield, - r_cut=r_cut, - dt=dt, - device=device, - seed=seed, - gsd_write_freq=gsd_write_freq, - gsd_file_name=gsd_file_name, - log_write_freq=log_write_freq, + initial_state=initial_state, + forcefield=forcefield, + r_cut=r_cut, + dt=dt, + device=device, + seed=seed, + gsd_write_freq=gsd_write_freq, + gsd_file_name=gsd_file_name, + log_write_freq=log_write_freq, + log_file_name=log_file_name, ) - self.tensile_axis=tensile_axis.lower() + self.tensile_axis = tensile_axis.lower() self.fix_ratio = fix_ratio axis_array_dict = { - "x": np.array([1,0,0]), - "y": np.array([0,1,0]), - "z": np.array([0,0,1]) + "x": np.array([1, 0, 0]), + "y": np.array([0, 1, 0]), + "z": np.array([0, 0, 1]), } axis_dict = {"x": 0, "y": 1, "z": 2} self._axis_index = axis_dict[self.tensile_axis] @@ -48,7 +49,7 @@ def __init__( self.fix_length = self.initial_length * fix_ratio # Set up walls of fixed particles: snapshot = self.state.get_snapshot() - positions = snapshot.particles.position[:,self._axis_index] + positions = snapshot.particles.position[:, self._axis_index] box_max = self.initial_length / 2 box_min = -box_max left_tags = np.where(positions < (box_min + self.fix_length))[0] @@ -58,44 +59,46 @@ def __init__( all_fixed = hoomd.filter.Union(self.fix_left, self.fix_right) # Set the group of particles to be integrated over self.integrate_group = hoomd.filter.SetDifference( - hoomd.filter.All(), all_fixed + hoomd.filter.All(), all_fixed ) @property def strain(self): - delta_L = self.box_lengths_reduced[self._axis_index]-self.initial_length + delta_L = ( + self.box_lengths_reduced[self._axis_index] - self.initial_length + ) return delta_L / self.initial_length def run_tensile(self, strain, kT, n_steps, period): current_length = self.box_lengths_reduced[self._axis_index] - final_length = current_length * (1+strain) + final_length = current_length * (1 + strain) final_box = np.copy(self.box_lengths_reduced) final_box[self._axis_index] = final_length - shift_by = (final_length - current_length) / (n_steps//period) + shift_by = (final_length - current_length) / (n_steps // period) resize_trigger = hoomd.trigger.Periodic(period) box_ramp = hoomd.variant.Ramp( - A=0, B=1, t_start=self.timestep, t_ramp=int(n_steps) + A=0, B=1, t_start=self.timestep, t_ramp=int(n_steps) ) box_resizer = hoomd.update.BoxResize( - box1=self.box_lengths_reduced, - box2=final_box, - variant=box_ramp, - trigger=resize_trigger, - filter=hoomd.filter.Null() + box1=self.box_lengths_reduced, + box2=final_box, + variant=box_ramp, + trigger=resize_trigger, + filter=hoomd.filter.Null(), ) particle_puller = PullParticles( - shift_by=shift_by/2, + shift_by=shift_by / 2, axis=self._axis_array, neg_filter=self.fix_left, - pos_filter=self.fix_right + pos_filter=self.fix_right, ) particle_updater = hoomd.update.CustomUpdater( - trigger=resize_trigger, action=particle_puller + trigger=resize_trigger, action=particle_puller ) self.operations.updaters.append(box_resizer) self.operations.updaters.append(particle_updater) self.set_integrator_method( integrator_method=hoomd.md.methods.NVE, - method_kwargs={"filter": self.integrate_group} + method_kwargs={"filter": self.integrate_group}, ) self.run(n_steps + 1) diff --git a/hoomd_polymers/library/systems.py b/hoomd_polymers/library/systems.py new file mode 100644 index 00000000..68ad997d --- /dev/null +++ b/hoomd_polymers/library/systems.py @@ -0,0 +1 @@ +# This is a placeholder for any class that inherits from base.system diff --git a/hoomd_polymers/modules/__init__.py b/hoomd_polymers/modules/__init__.py new file mode 100644 index 00000000..9ab8ed42 --- /dev/null +++ b/hoomd_polymers/modules/__init__.py @@ -0,0 +1 @@ +from .welding import Interface, SlabSimulation, WeldSimulation diff --git a/hoomd_polymers/modules/utils.py b/hoomd_polymers/modules/utils.py index fffb2d20..7fe5efa1 100644 --- a/hoomd_polymers/modules/utils.py +++ b/hoomd_polymers/modules/utils.py @@ -1,21 +1,14 @@ -import gsd.hoomd import hoomd import numpy as np def add_void_particles( - snapshot, - forcefield, - num_voids, - void_axis, - void_diameter, - epsilon, - r_cut - ): + snapshot, forcefield, num_voids, void_axis, void_diameter, epsilon, r_cut +): void_axis = np.asarray(void_axis) snapshot.particles.N += num_voids snapshot.particles.position[-1] = ( - void_axis * snapshot.configuration.box[0:3]/2 + void_axis * snapshot.configuration.box[0:3] / 2 ) snapshot.particles.types = snapshot.particles.types + ["VOID"] snapshot.particles.typeid[-1] = len(snapshot.particles.types) - 1 @@ -23,7 +16,8 @@ def add_void_particles( lj = [i for i in forcefield if isinstance(i, hoomd.md.pair.LJ)][0] for ptype in snapshot.particles.types: lj.params[(ptype, "VOID")] = { - 'sigma': void_diameter, 'epsilon': epsilon + "sigma": void_diameter, + "epsilon": epsilon, } lj.r_cut[(ptype, "VOID")] = r_cut diff --git a/hoomd_polymers/modules/welding.py b/hoomd_polymers/modules/welding.py index 30e4a3d8..7b9f3ff7 100644 --- a/hoomd_polymers/modules/welding.py +++ b/hoomd_polymers/modules/welding.py @@ -2,8 +2,7 @@ import hoomd import numpy as np -from hoomd_polymers.sim.simulation import Simulation -from hoomd_polymers.systems import System +from hoomd_polymers.base.simulation import Simulation class Interface: @@ -31,63 +30,72 @@ def _build(self): interface.dihedrals.M = snap.dihedrals.M interface.pairs.N = snap.pairs.N * 2 - # Set up box. Box edge is doubled along the interface axis direction, plus the gap + # Set up box. Box edge is doubled along the interface axis direction, + # plus the gap interface.configuration.box = np.copy(snap.configuration.box) interface.configuration.box[axis_index] *= 2 - interface.configuration.box[axis_index] += (self.gap - self.wall_sigma) - + interface.configuration.box[axis_index] += self.gap - self.wall_sigma + # Set up snapshot.particles info: # Get set of new coordiantes, shifted along interface axis - shift = (snap.configuration.box[axis_index]+self.gap-self.wall_sigma)/2 + shift = ( + snap.configuration.box[axis_index] + self.gap - self.wall_sigma + ) / 2 right_pos = np.copy(snap.particles.position) - right_pos[:,axis_index] += shift + right_pos[:, axis_index] += shift left_pos = np.copy(snap.particles.position) - left_pos[:,axis_index] -= shift - + left_pos[:, axis_index] -= shift + pos = np.concatenate((left_pos, right_pos), axis=None) - mass = np.concatenate((snap.particles.mass, snap.particles.mass), axis=None) - charges = np.concatenate((snap.particles.charge, snap.particles.charge), axis=None) + mass = np.concatenate( + (snap.particles.mass, snap.particles.mass), axis=None + ) + charges = np.concatenate( + (snap.particles.charge, snap.particles.charge), axis=None + ) type_ids = np.concatenate( - (snap.particles.typeid, snap.particles.typeid), axis=None + (snap.particles.typeid, snap.particles.typeid), axis=None ) interface.particles.position = pos interface.particles.mass = mass interface.particles.charge = charges interface.particles.types = snap.particles.types interface.particles.typeid = type_ids - + # Set up bonds: bond_group_left = np.copy(snap.bonds.group) bond_group_right = np.copy(snap.bonds.group) + snap.particles.N - bond_group = np.concatenate((bond_group_left, bond_group_right), axis=None) + bond_group = np.concatenate( + (bond_group_left, bond_group_right), axis=None + ) bond_type_ids = np.concatenate( - (snap.bonds.typeid, snap.bonds.typeid), axis=None + (snap.bonds.typeid, snap.bonds.typeid), axis=None ) interface.bonds.group = bond_group interface.bonds.typeid = bond_type_ids interface.bonds.types = snap.bonds.types - + # Set up angles: angle_group_left = np.copy(snap.angles.group) angle_group_right = np.copy(snap.angles.group) + snap.particles.N angle_group = np.concatenate( - (angle_group_left, angle_group_right), axis=None + (angle_group_left, angle_group_right), axis=None ) angle_type_ids = np.concatenate( - (snap.angles.typeid, snap.angles.typeid), axis=None + (snap.angles.typeid, snap.angles.typeid), axis=None ) interface.angles.group = angle_group interface.angles.typeid = angle_type_ids interface.angles.types = snap.angles.types - + # Set up dihedrals: dihedral_group_left = np.copy(snap.dihedrals.group) dihedral_group_right = np.copy(snap.dihedrals.group) + snap.particles.N dihedral_group = np.concatenate( - (dihedral_group_left, dihedral_group_right), axis=None + (dihedral_group_left, dihedral_group_right), axis=None ) dihedral_type_ids = np.concatenate( - (snap.dihedrals.typeid, snap.dihedrals.typeid), axis=None + (snap.dihedrals.typeid, snap.dihedrals.typeid), axis=None ) interface.dihedrals.group = dihedral_group interface.dihedrals.typeid = dihedral_type_ids @@ -99,7 +107,7 @@ def _build(self): pair_group_right = np.copy(snap.pairs.group) + snap.particles.N pair_group = np.concatenate((pair_group_left, pair_group_right)) pair_type_ids = np.concatenate( - (snap.pairs.typeid, snap.pairs.typeid), axis=None + (snap.pairs.typeid, snap.pairs.typeid), axis=None ) interface.pairs.group = pair_group interface.pairs.typeid = pair_type_ids @@ -109,42 +117,42 @@ def _build(self): class SlabSimulation(Simulation): def __init__( - self, - initial_state, - forcefield, - interface_axis="x", - wall_sigma=1.0, - wall_epsilon=1.0, - wall_r_cut=2.5, - wall_r_extrap=0, - r_cut=2.5, - seed=42, - gsd_write_freq=1e4, - gsd_file_name="weld.gsd", - log_write_freq=1e3, - log_file_name="sim_data.txt" - ): + self, + initial_state, + forcefield, + interface_axis="x", + wall_sigma=1.0, + wall_epsilon=1.0, + wall_r_cut=2.5, + wall_r_extrap=0, + r_cut=2.5, + seed=42, + gsd_write_freq=1e4, + gsd_file_name="weld.gsd", + log_write_freq=1e3, + log_file_name="sim_data.txt", + ): super(SlabSimulation, self).__init__( - initial_state=initial_state, - forcefield=forcefield, - r_cut=r_cut, - seed=seed, - gsd_write_freq=gsd_write_freq, - gsd_file_name=gsd_file_name, - log_write_freq=log_write_freq, - log_file_name=log_file_name + initial_state=initial_state, + forcefield=forcefield, + r_cut=r_cut, + seed=seed, + gsd_write_freq=gsd_write_freq, + gsd_file_name=gsd_file_name, + log_write_freq=log_write_freq, + log_file_name=log_file_name, ) self.interface_axis = interface_axis.lower() - axis_array_dict = {"x": (1,0,0), "y": (0, 1, 0), "z": (0, 0, 1)} + axis_array_dict = {"x": (1, 0, 0), "y": (0, 1, 0), "z": (0, 0, 1)} axis_dict = {"x": 0, "y": 1, "z": 2} self._axis_array = axis_array_dict[self.interface_axis] self._axis_index = axis_dict[self.interface_axis] self.add_walls( - self._axis_array, - wall_sigma, - wall_epsilon, - wall_r_cut, - wall_r_extrap + self._axis_array, + wall_sigma, + wall_epsilon, + wall_r_cut, + wall_r_extrap, ) snap = self.state.get_snapshot() @@ -154,38 +162,34 @@ def __init__( class WeldSimulation(Simulation): def __init__( - self, - initial_state, - forcefield, - interface_axis="x", - wall_sigma=1.0, - wall_epsilon=1.0, - wall_r_cut=2.5, - wall_r_extrap=0, - r_cut=2.5, - seed=42, - gsd_write_freq=1e4, - gsd_file_name="weld.gsd", - log_write_freq=1e3, - log_file_name="sim_data.txt" + self, + initial_state, + forcefield, + interface_axis="x", + wall_sigma=1.0, + wall_epsilon=1.0, + wall_r_cut=2.5, + wall_r_extrap=0, + r_cut=2.5, + seed=42, + gsd_write_freq=1e4, + gsd_file_name="weld.gsd", + log_write_freq=1e3, + log_file_name="sim_data.txt", ): super(WeldSimulation, self).__init__( - initial_state=initial_state, - forcefield=forcefield, - r_cut=r_cut, - seed=seed, - gsd_write_freq=gsd_write_freq, - gsd_file_name=gsd_file_name, - log_write_freq=log_write_freq, - log_file_name=log_file_name + initial_state=initial_state, + forcefield=forcefield, + r_cut=r_cut, + seed=seed, + gsd_write_freq=gsd_write_freq, + gsd_file_name=gsd_file_name, + log_write_freq=log_write_freq, + log_file_name=log_file_name, ) - axis_dict = {"x": (1,0,0), "y": (0, 1, 0), "z": (0, 0, 1)} + axis_dict = {"x": (1, 0, 0), "y": (0, 1, 0), "z": (0, 0, 1)} self.interface_axis = interface_axis.lower() self.wall_axis = axis_dict[self.interface_axis] self.add_walls( - self.wall_axis, - wall_sigma, - wall_epsilon, - wall_r_cut, - wall_r_extrap + self.wall_axis, wall_sigma, wall_epsilon, wall_r_cut, wall_r_extrap ) diff --git a/hoomd_polymers/molecules.py b/hoomd_polymers/molecules.py deleted file mode 100644 index 875fc496..00000000 --- a/hoomd_polymers/molecules.py +++ /dev/null @@ -1,253 +0,0 @@ -import os -import random - -import mbuild as mb -from mbuild.coordinate_transform import z_axis_transform -from mbuild.lib.recipes import Polymer -import numpy as np - -from hoomd_polymers.library import MON_DIR - - -class CoPolymer(Polymer): - """Builds a polymer consisting of two monomer types. - - Parameters - ---------- - monomer_A : hoomd_polymers.molecules.Polymer; required - Class of the A-type monomer - monomer_B : hoomd_polymers.molecules.Polymer: required - Class of the B-type monomer - length : int; required - The total number of monomers in the molecule - sequence : str; optional; default None - Manually define the sequence of 'A' and 'B' monomers. - Leave as None if generating random sequences. - Example: sequence = "AABAABAAB" - random_sequence : bool; optional; default True - Creates a random 'A' 'B' sequence as a function - of the AB_ratio. Set to False when manually - defining sequence - AB_ratio : float; optional; default 0.50 - The relative weight of A to B monomer types. - Used when generating random sequences. - seed : int; optional; default 24 - Set the seed used when generating random sequences - """ - def __init__( - self, - monomer_A, - monomer_B, - length, - sequence=None, - random_sequence=True, - AB_ratio=0.50, - seed=24 - ): - super(CoPolymer, self).__init__() - self.monomer_A = monomer_A(length=1) - self.monomer_B = monomer_B(length=1) - if random_sequence: - random.seed(seed) - self.sequence = random.choices( - ["A", "B"], [AB_ratio, 1-AB_ratio], k=length - ) - length = 1 - else: - self.sequence = sequence - self.A_ratio = self.sequence.count("A")/len(self.sequence) - self.B_ratio = self.sequence.count("B")/len(self.sequence) - - self.add_monomer( - self.monomer_A.monomer, - indices=self.monomer_A.bond_indices, - orientation=self.monomer_A.bond_orientation, - separation=self.monomer_A.bond_length - ) - self.add_monomer( - self.monomer_B.monomer, - indices=self.monomer_B.bond_indices, - orientation=self.monomer_B.bond_orientation, - separation=self.monomer_B.bond_length - ) - self.build(n=length, sequence=self.sequence) - - -class PolyEthylene(Polymer): - """Creates a Poly(ethylene) chain. - - Parameters - ---------- - length : int; required - The number of monomer repeat units in the chain - """ - def __init__(self, length): - super(PolyEthylene, self).__init__() - self.smiles_str = "CC" - self.description = "Poly(ethylene)" - self.file = None - self.monomer = mb.load(self.smiles_str, smiles=True) - self.bond_indices = [2, 6] - self.bond_length = 0.145 - self.bond_orientation = [None, None] - self.add_monomer( - self.monomer, - indices=self.bond_indices, - separation=self.bond_length - ) - self.build(n=length, sequence="A") - # Align the chain along the z-axis - z_axis_transform( - self, - point_on_z_axis=self[-2], - point_on_zx_plane=self[-1] - ) - - -class PPS(Polymer): - """Creates a Poly(phenylene-sulfide) (PPS) chain. - - Parameters - ---------- - length : int; required - The number of monomer repeat units in the chain - """ - def __init__(self, length): - super(PPS, self).__init__() - self.smiles_str = "c1ccc(S)cc1" - self.file = None - self.description = "Poly(phenylene-sulfide)" - self.monomer = mb.load(self.smiles_str, smiles=True) - # Need to align monomer along zx plane due to orientation of S-H bond - z_axis_transform( - self.monomer, - point_on_z_axis=self.monomer[7], - point_on_zx_plane=self.monomer[4] - ) - self.bond_indices = [7, 10] - self.bond_length = 0.176 - self.bond_orientation = [[0, 0, 1], [0, 0, -1]] - self.add_monomer( - self.monomer, - indices=self.bond_indices, - separation=self.bond_length, - orientation=self.bond_orientation - ) - self.build(n=length, sequence="A") - - -class PEEK(Polymer): - def __init__(self, length): - super(PEEK, self).__init__() - - -class PEKK_para(Polymer): - """Creates a Poly(ether-ketone-ketone) (PEKK) chain. - The bonding positions of consecutive ketone groups - takes place on the para site of the phenyl ring. - - Parameters - ---------- - length : int; required - The number of monomer repeat units in the chain - """ - def __init__(self, length): - super(PEKK_para, self).__init__() - self.smiles_str = "c1ccc(Oc2ccc(C(=O)c3ccc(C(=O))cc3)cc2)cc1" - self.file = os.path.join(MON_DIR, "pekk_para.mol2") - self.description = ("Poly(ether-ketone-ketone) with para bonding " - "configuration between consectuvie " - "ketone linkage groups") - self.monomer = mb.load(self.file) - self.bond_indices = [35, 36] - self.bond_length = 0.148 - self.bond_orientation = [[0, 0, -1], [0, 0, 1]] - self.add_monomer( - self.monomer, - indices=self.bond_indices, - separation=self.bond_length, - orientation=self.bond_orientation - ) - self.build(n=length, sequence="A") - - -class PEKK_meta(Polymer): - """Creates a Poly(ether-ketone-ketone) (PEKK) chain. - The bonding positions of consecutive ketone groups - takes place on the meta site of the phenyl ring. - - Parameters - ---------- - length : int; required - The number of monomer repeat units in the chain - """ - def __init__(self, length): - super(PEKK_meta, self).__init__() - self.smiles_str = "c1cc(Oc2ccc(C(=O)c3cc(C(=O))ccc3)cc2)ccc1" - self.file = os.path.join(MON_DIR, "pekk_meta.mol2") - self.description = ("Poly(ether-ketone-ketone) with meta bonding " - "configuration between consectuvie " - "ketone linkage groups") - self.monomer = mb.load(self.file) - self.bond_indices = [35, 36] - self.bond_length = 0.148 - self.bond_orientation = [[0, 0, -1], [0, 0, 1]] - self.add_monomer( - self.monomer, - indices=self.bond_indices, - separation=self.bond_length, - orientation=self.bond_orientation - ) - self.build(n=length, sequence="A") - - -class LJ_chain(mb.Compound): - """Creates a coarse-grained bead-spring polymer chain. - - Parameters - ---------- - length : int; required - The number of times to repeat bead_sequence in a single chain. - bead_sequence : list; optional; default ["A"] - The sequence of bead types in the chain. - bond_length : dict; optional; default {"A-A": 1.0} - The bond length between connected beads (units: nm) - bead_mass : dict; optional; default {"A": 1.0} - The mass of the bead types - """ - def __init__( - self, - length, - bead_sequence=["A"], - bead_mass={"A": 1.0}, - bond_lengths={"A-A": 1.0}, - ): - super(LJ_chain, self).__init__() - self.description = "Simple bead-spring polymer" - last_bead = None - for i in range(length): - for idx, bead_type in enumerate(bead_sequence): - mass = bead_mass.get(bead_type, None) - if not mass: - raise ValueError( - f"The bead mass for {bead_type} was not given " - "in the bead_mass dict." - ) - next_bead = mb.Compound(mass=mass, name=bead_type, charge=0) - self.add(next_bead) - if last_bead: - bead_pair = "-".join([last_bead.name, next_bead.name]) - bond_length = bond_lengths.get(bead_pair, None) - if not bond_length: - bead_pair_rev = "-".join([next_bead.name, last_bead.name]) - bond_length = bond_lengths.get(bead_pair_rev, None) - if not bond_length: - raise ValueError( - "The bond length for pair " - f"{bead_pair} or {bead_pair_rev} " - "is not found in the bond_lengths dict." - ) - new_pos = last_bead.xyz[0] + (0, 0, bond_length) - next_bead.translate_to(new_pos) - self.add_bond([next_bead, last_bead]) - last_bead = next_bead diff --git a/hoomd_polymers/sim/__init__.py b/hoomd_polymers/sim/__init__.py deleted file mode 100644 index eb29d0aa..00000000 --- a/hoomd_polymers/sim/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -from .simulation import Simulation -from .tensile import Tensile diff --git a/hoomd_polymers/systems.py b/hoomd_polymers/systems.py deleted file mode 100644 index eac10a2e..00000000 --- a/hoomd_polymers/systems.py +++ /dev/null @@ -1,312 +0,0 @@ -import mbuild as mb -from mbuild.formats.hoomd_forcefield import create_hoomd_forcefield -import numpy as np -import unyt - - -from hoomd_polymers.utils import scale_charges, check_return_iterable - - -class System: - """Base class from which other systems inherit. - - Parameters - ---------- - molecule : hoomd_polymers.molecule; required - n_mols : int; required - The number of times to replicate molecule in the system - density : float; optional; default None - The desired density of the system (g/cm^3). Used to set the - target_box attribute. Can be useful when initializing - systems at low denisty and running a shrink simulaton - to acheive a target density. - """ - def __init__(self, molecule, n_mols, mol_kwargs={}, density=None): - self.n_mols = check_return_iterable(n_mols) - self._molecules = check_return_iterable(molecule) - self.mol_kwargs = check_return_iterable(mol_kwargs) - self.density = density - self.target_box = None - self.system = None - self.typed_system = None - self._hoomd_objects = None - self._reference_values = None - self.molecules = [] - - for mol, n, kw_args, in zip( - self._molecules, - self.n_mols, - self.mol_kwargs - ): - for i in range(n): - self.molecules.append(mol(**kw_args)) - - @property - def mass(self): - if not self.system: - return sum(i.mass for i in self.molecules) - else: - return self.system.mass - - @property - def box(self): - return self.system.box - - @property - def hoomd_snapshot(self): - if not self._hoomd_objects: - raise ValueError( - "The hoomd snapshot has not yet been created. " - "Create a Hoomd snapshot and forcefield by applying " - "a forcefield using System.apply_forcefield()." - ) - else: - return self._hoomd_objects[0] - - @property - def hoomd_forcefield(self): - if not self._hoomd_objects: - raise ValueError( - "The hoomd forcefield has not yet been created. " - "Create a Hoomd snapshot and forcefield by applying " - "a forcefield using System.apply_forcefield()." - ) - else: - return self._hoomd_objects[1] - - @property - def reference_distance(self): - return self._reference_values.distance * unyt.angstrom - - @property - def reference_mass(self): - return self._reference_values.mass * unyt.amu - - @property - def reference_energy(self): - return self._reference_values.energy * unyt.kcal / unyt.mol - - def apply_forcefield( - self, - forcefield, - remove_hydrogens=False, - scale_parameters=True, - remove_charges=False, - make_charge_neutral=False, - r_cut=2.5 - ): - if len(self._molecules) == 1: - use_residue_map = True - else: - use_residue_map = False - self.typed_system = forcefield.apply( - structure=self.system, use_residue_map=use_residue_map - ) - if remove_hydrogens: - print("Removing hydrogen atoms and adjusting heavy atoms") - # Try by element first: - hydrogens = [a for a in self.typed_system.atoms if a.element == 1] - if len(hydrogens) == 0: # Try by mass - hydrogens = [a for a in self.typed_system.atoms if a.mass == 1.008] - if len(hydrogens) == 0: - warnings.warn( - "Hydrogen atoms could not be found by element or mass" - ) - for h in hydrogens: - h.atomic_number = 1 - bonded_atom = h.bond_partners[0] - bonded_atom.mass += h.mass - bonded_atom.charge += h.charge - self.typed_system.strip( - [a.atomic_number == 1 for a in self.typed_system.atoms] - ) - if remove_charges: - for atom in self.typed_system.atoms: - atom.charge = 0 - if make_charge_neutral and not remove_charges: - print("Adjust charges to make system charge neutral") - new_charges = scale_charges( - charges=np.array([a.charge for a in self.typed_system.atoms]), - n_particles=len(self.typed_system.atoms) - ) - for idx, charge in enumerate(new_charges): - self.typed_system.atoms[idx].charge = charge - - init_snap, forcefield, refs = create_hoomd_forcefield( - structure=self.typed_system, - r_cut=r_cut, - auto_scale=scale_parameters - ) - self._hoomd_objects = [init_snap, forcefield] - self._reference_values = refs - - def set_target_box( - self, x_constraint=None, y_constraint=None, z_constraint=None - ): - """Set the target volume of the system during - the initial shrink step. - If no constraints are set, the target box is cubic. - Setting constraints will hold those box vectors - constant and adjust others to match the target density. - - Parameters - ----------- - x_constraint : float, optional, defualt=None - Fixes the box length (nm) along the x axis - y_constraint : float, optional, default=None - Fixes the box length (nm) along the y axis - z_constraint : float, optional, default=None - Fixes the box length (nm) along the z axis - - """ - if not any([x_constraint, y_constraint, z_constraint]): - Lx = Ly = Lz = self._calculate_L() - else: - constraints = np.array([x_constraint, y_constraint, z_constraint]) - fixed_L = constraints[np.where(constraints!=None)] - #Conv from nm to cm for _calculate_L - fixed_L *= 1e-7 - L = self._calculate_L(fixed_L = fixed_L) - constraints[np.where(constraints==None)] = L - Lx, Ly, Lz = constraints - - self.target_box = np.array([Lx, Ly, Lz]) - - def visualize(self): - if self.system: - self.system.visualize().show() - else: - raise ValueError( - "The initial configuraiton has not been created yet." - ) - - def _calculate_L(self, fixed_L=None): - """Calculates the required box length(s) given the - mass of a sytem and the target density. - - Box edge length constraints can be set by set_target_box(). - If constraints are set, this will solve for the required - lengths of the remaining non-constrained edges to match - the target density. - - Parameters - ---------- - fixed_L : np.array, optional, defualt=None - Array of fixed box lengths to be accounted for - when solving for L - - """ - # Convert from amu to grams - M = self.mass * 1.66054e-24 - vol = (M / self.density) # cm^3 - if fixed_L is None: - L = vol**(1/3) - else: - L = vol / np.prod(fixed_L) - if len(fixed_L) == 1: # L is cm^2 - L = L**(1/2) - # Convert from cm back to nm - L *= 1e7 - return L - - -class Pack(System): - """Uses PACKMOL via mbuild.packing.fill_box. - The box used for packing is expanded to allow PACKMOL to place all of the molecules. - - Parameters - ---------- - packing_expand_factor : int; optional, default 5 - - """ - def __init__( - self, - molecule, - n_mols, - mol_kwargs={}, - density=None, - packing_expand_factor=5, - edge=0.2 - ): - super(Pack, self).__init__( - molecule=molecule, - n_mols=n_mols, - mol_kwargs=mol_kwargs, - density=density - ) - self.packing_expand_factor = packing_expand_factor - self.edge = edge - self._build() - - def _build(self): - self.set_target_box() - self.system = mb.packing.fill_box( - compound=self.molecules, - n_compounds=[1 for i in self.molecules], - box=list(self.target_box*self.packing_expand_factor), - overlap=0.2, - edge=self.edge - ) - - -class Lattice(System): - """Places the molecules in a lattice configuration. - Assumes two molecules per unit cell. - - Parameters - ---------- - x : float; required - The distance (nm) between lattice points in the x direction. - y : float; required - The distance (nm) between lattice points in the y direction. - n : int; required - The number of times to repeat the unit cell in x and y - lattice_vector : array-like - The vector between points in the unit cell - """ - def __init__( - self, - molecule, - density, - n_mols, - x, - y, - n, - mol_kwargs={}, - basis_vector=[0.5, 0.5, 0], - z_adjust=1.0, - ): - super(Lattice, self).__init__( - molecule=molecule, - n_mols=n_mols, - mol_kwargs=mol_kwargs, - density=density - ) - self.x = x - self.y = y - self.n = n - self.basis_vector = basis_vector - self._build() - - def _build(self): - next_idx = 0 - self.system = mb.Compound() - for i in range(self.n): - layer = mb.Compound() - for j in range(self.n): - try: - comp1 = self.molecules[next_idx] - comp2 = self.molecules[next_idx + 1] - comp2.translate(self.basis_vector) - unit_cell = mb.Compound(subcompounds=[comp1, comp2]) - unit_cell.translate((0, self.y*j, 0)) - layer.add(unit_cell) - next_idx += 2 - except IndexError: - pass - layer.translate((self.x*i, 0, 0)) - self.system.add(layer) - bounding_box = self.system.get_boundingbox() - x_len = bounding_box.lengths[0] - y_len = bounding_box.lengths[1] - self.set_target_box(x_constraint=x_len, y_constraint=y_len) diff --git a/hoomd_polymers/tests/__init__.py b/hoomd_polymers/tests/__init__.py new file mode 100644 index 00000000..9eb4c5fb --- /dev/null +++ b/hoomd_polymers/tests/__init__.py @@ -0,0 +1 @@ +from .base_test import BaseTest diff --git a/hoomd_polymers/tests/assets/benzene.mol2 b/hoomd_polymers/tests/assets/benzene.mol2 new file mode 100644 index 00000000..5f30f1ad --- /dev/null +++ b/hoomd_polymers/tests/assets/benzene.mol2 @@ -0,0 +1,32 @@ +@MOLECULE +***** + 12 12 0 0 0 +SMALL +GASTEIGER + +@ATOM + 1 C -0.1172 2.3334 -0.0000 C.ar 1 UNL1 -0.0618 + 2 C -1.3248 1.6367 -0.0000 C.ar 1 UNL1 -0.0618 + 3 C -1.3252 0.2425 -0.0000 C.ar 1 UNL1 -0.0618 + 4 C -0.1180 -0.4549 -0.0000 C.ar 1 UNL1 -0.0618 + 5 C 1.0896 0.2418 -0.0000 C.ar 1 UNL1 -0.0618 + 6 C 1.0900 1.6360 -0.0000 C.ar 1 UNL1 -0.0618 + 7 H -0.1169 3.4207 -0.0000 H 1 UNL1 0.0618 + 8 H -2.2663 2.1806 -0.0000 H 1 UNL1 0.0618 + 9 H -2.2670 -0.3009 -0.0000 H 1 UNL1 0.0618 + 10 H -0.1183 -1.5422 -0.0000 H 1 UNL1 0.0618 + 11 H 2.0310 -0.3021 -0.0000 H 1 UNL1 0.0618 + 12 H 2.0318 2.1793 -0.0000 H 1 UNL1 0.0618 +@BOND + 1 1 2 ar + 2 2 3 ar + 3 3 4 ar + 4 4 5 ar + 5 5 6 ar + 6 1 6 ar + 7 1 7 1 + 8 2 8 1 + 9 3 9 1 + 10 4 10 1 + 11 5 11 1 + 12 6 12 1 diff --git a/hoomd_polymers/tests/assets/benzene_oplsaa.xml b/hoomd_polymers/tests/assets/benzene_oplsaa.xml new file mode 100644 index 00000000..437f644a --- /dev/null +++ b/hoomd_polymers/tests/assets/benzene_oplsaa.xml @@ -0,0 +1,23 @@ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/hoomd_polymers/tests/assets/test_ff.xml b/hoomd_polymers/tests/assets/test_ff.xml new file mode 100644 index 00000000..ccaa19e7 --- /dev/null +++ b/hoomd_polymers/tests/assets/test_ff.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/hoomd_polymers/tests/base/test_molecule.py b/hoomd_polymers/tests/base/test_molecule.py new file mode 100644 index 00000000..53cbb5f9 --- /dev/null +++ b/hoomd_polymers/tests/base/test_molecule.py @@ -0,0 +1,284 @@ +import pytest + +from hoomd_polymers import CoPolymer, Molecule, Polymer +from hoomd_polymers.tests import BaseTest +from hoomd_polymers.utils import FF_Types, exceptions + + +class TestMolecule(BaseTest): + def test_molecule_from_mb_compound(self, benzene_mb): + molecule = Molecule(num_mols=2, compound=benzene_mb) + assert len(molecule.molecules) == 2 + + def test_molecule_from_gmso_topology(self, benzene_gmso): + molecule = Molecule(num_mols=2, compound=benzene_gmso) + assert len(molecule.molecules) == 2 + + def test_molecule_from_smiles(self, benzene_smiles): + molecule = Molecule(num_mols=2, smiles=benzene_smiles) + assert len(molecule.molecules) == 2 + + def test_molecule_from_file(self, benzene_mol2): + molecule = Molecule(num_mols=2, file=benzene_mol2) + assert len(molecule.molecules) == 2 + + def test_molecule_topology_benzene(self, benzene_mb): + molecule = Molecule(num_mols=2, compound=benzene_mb) + assert set(molecule.topology_information["particle_types"]) == { + "C", + "H", + } + assert set(molecule.topology_information["pair_types"]) == { + ("C", "C"), + ("C", "H"), + ("H", "H"), + } + assert len(set(molecule.topology_information["particle_typeid"])) == 2 + assert len(molecule.topology_information["bond_types"]) == 2 + assert len(molecule.topology_information["angle_types"]) == 2 + assert len(molecule.topology_information["dihedral_types"]) == 3 + assert not any(molecule.topology_information["particle_charge"]) + + def test_validate_force_field_oplsaa(self, benzene_mb): + molecule = Molecule( + num_mols=2, force_field="oplsaa", compound=benzene_mb + ) + assert molecule.ff_type == FF_Types.oplsaa + assert set(molecule.topology_information["particle_types"]) == { + "opls_145", + "opls_146", + } + assert any(molecule.topology_information["particle_charge"]) + + def test_validate_force_field_xml_file(self, benzene_mb): + molecule = Molecule( + num_mols=2, force_field="oplsaa.xml", compound=benzene_mb + ) + assert molecule.ff_type == FF_Types.oplsaa + assert set(molecule.topology_information["particle_types"]) == { + "opls_145", + "opls_146", + } + assert any(molecule.topology_information["particle_charge"]) + + def test_validate_force_field_xml_file_path(self, benzene_mb, benzene_xml): + molecule = Molecule( + num_mols=2, force_field=benzene_xml, compound=benzene_mb + ) + assert molecule.ff_type == FF_Types.custom + assert set(molecule.topology_information["particle_types"]) == { + "opls_145", + "opls_146", + } + assert any(molecule.topology_information["particle_charge"]) + + def test_validate_force_field_not_xml_file(self, benzene_mb): + with pytest.raises(ValueError): + Molecule(num_mols=2, force_field="oplsaa.txt", compound=benzene_mb) + + def test_validate_force_field_not_supported(self, benzene_mb): + with pytest.raises(ValueError): + Molecule(num_mols=2, force_field="oplsaa2", compound=benzene_mb) + + def test_validate_force_field_invalid_xml_file(self, benzene_mb): + with pytest.raises(ValueError): + Molecule(num_mols=2, force_field="oplsaa2.xml", compound=benzene_mb) + + def test_validate_force_field_hoomd_ff_aa( + self, benzene_mb, benzene_hoomd_ff + ): + hoomd_ff = benzene_hoomd_ff(include_hydrogen=True) + molecule = Molecule( + num_mols=2, force_field=hoomd_ff, compound=benzene_mb + ) + assert molecule.ff_type == FF_Types.Hoomd + + def test_validate_fore_field_hoomd_ff_ua( + self, benzene_mb, benzene_hoomd_ff + ): + hoomd_ff = benzene_hoomd_ff(include_hydrogen=False) + molecule = Molecule( + num_mols=2, force_field=hoomd_ff, compound=benzene_mb + ) + assert molecule.ff_type == FF_Types.Hoomd + + def test_validate_force_field_hoomd_ff_missing_pair( + self, benzene_mb, benzene_hoomd_ff + ): + hoomd_ff = benzene_hoomd_ff(include_hydrogen=True) + hoomd_ff.pop(0) + with pytest.raises(exceptions.MissingPairPotentialError): + Molecule(num_mols=2, force_field=hoomd_ff, compound=benzene_mb) + + def test_validate_force_field_hoomd_ff_missing_bond( + self, benzene_mb, benzene_hoomd_ff + ): + hoomd_ff = benzene_hoomd_ff(include_hydrogen=True) + hoomd_ff.pop(1) + with pytest.raises(exceptions.MissingBondPotentialError): + Molecule(num_mols=2, force_field=hoomd_ff, compound=benzene_mb) + + def test_validate_force_field_hoomd_ff_invalid_pair( + self, benzene_mb, benzene_hoomd_ff + ): + hoomd_ff = benzene_hoomd_ff(include_hydrogen=True, invalid_pair=True) + with pytest.raises(exceptions.MissingPairPotentialError): + Molecule(num_mols=2, force_field=hoomd_ff, compound=benzene_mb) + + def test_validate_force_field_hoomd_ff_missing_Coulomb( + self, benzene_mb, benzene_xml, benzene_hoomd_ff + ): + hoomd_ff = benzene_hoomd_ff(include_hydrogen=True) + typed_molecule = Molecule( + num_mols=2, force_field=benzene_xml, compound=benzene_mb + ) + with pytest.raises(exceptions.MissingCoulombPotentialError): + Molecule( + num_mols=2, + force_field=hoomd_ff, + compound=typed_molecule.gmso_molecule, + ) + + def test_coarse_grain_with_single_beads(self, benzene_smiles): + molecule = Molecule(num_mols=2, smiles=benzene_smiles) + molecule.coarse_grain(beads={"A": benzene_smiles}) + assert molecule.topology_information["particle_types"] == ["A"] + assert molecule.n_particles == 2 + assert molecule.n_bonds == 0 + + def test_coarse_grain_with_multiple_single_beads( + self, octane_smiles, ethane_smiles + ): + molecule = Molecule(num_mols=2, smiles=octane_smiles) + molecule.coarse_grain(beads={"A": ethane_smiles}) + assert molecule.molecules[0].n_particles == 4 + assert molecule.topology_information["bond_types"] == {("A", "A")} + + def test_coarse_grain_with_different_beads( + self, pps_smiles, benzene_smiles + ): + molecule = Molecule(num_mols=2, smiles=pps_smiles) + molecule.coarse_grain(beads={"A": benzene_smiles, "B": "S"}) + assert molecule.molecules[0].n_particles == 2 + assert molecule.topology_information["particle_types"] == ["A", "B"] + assert molecule.topology_information["bond_types"] == {("A", "B")} + + def test_coarse_grain_invalid_beads(self, benzene_smiles): + molecule = Molecule(num_mols=2, smiles=benzene_smiles) + with pytest.raises(ValueError): + molecule.coarse_grain(beads={"A": "CO"}) + + +class TestPolymer(BaseTest): + def test_polymer(self, dimethylether_smiles): + polymer = Polymer( + lengths=3, + num_mols=1, + smiles=dimethylether_smiles, + bond_indices=[3, -1], + bond_length=0.15, + bond_orientation=[None, None], + ) + assert polymer.n_particles == 23 + assert polymer.n_bonds == 22 + assert ("O", "C", "C") in polymer.topology_information["angle_types"] + assert ("O", "C", "C", "O") in polymer.topology_information[ + "dihedral_types" + ] + + def test_polymer_different_chain_lengths(self, dimethylether_smiles): + polymer = Polymer( + lengths=[3, 4], + num_mols=[1, 1], + smiles=dimethylether_smiles, + bond_indices=[3, -1], + bond_length=0.15, + bond_orientation=[None, None], + ) + assert polymer.n_particles == 53 + assert len(polymer.molecules[0].labels["monomer"]) == 3 + assert len(polymer.molecules[1].labels["monomer"]) == 4 + + def test_polymer_different_num_mol(self, dimethylether_smiles): + polymer = Polymer( + lengths=[3, 2], + num_mols=[1, 2], + smiles=dimethylether_smiles, + bond_indices=[3, -1], + bond_length=0.15, + bond_orientation=[None, None], + ) + assert polymer.n_particles == 55 + assert len(polymer.molecules[0].labels["monomer"]) == 3 + assert len(polymer.molecules[1].labels["monomer"]) == 2 + assert len(polymer.molecules[2].labels["monomer"]) == 2 + + def test_polymer_unequal_num_mol_length(self, dimethylether_smiles): + with pytest.raises(ValueError): + Polymer( + lengths=[3], + num_mols=[1, 2], + smiles=dimethylether_smiles, + bond_indices=[3, -1], + bond_length=0.15, + bond_orientation=[None, None], + ) + + +class TestCopolymer(BaseTest): + def test_copolymer_with_sequence(self, polyethylene, polyDME): + copolymer = CoPolymer( + monomer_A=polyDME, + monomer_B=polyethylene, + lengths=1, + num_mols=1, + sequence="ABA", + ) + assert copolymer.n_particles == 22 + assert ("C", "C", "C", "C") in copolymer.topology_information[ + "dihedral_types" + ] + + def test_copolymer_with_sequence_different_chain_lengths( + self, polyethylene, polyDME + ): + copolymer = CoPolymer( + monomer_A=polyDME, + monomer_B=polyethylene, + lengths=[2, 3], + num_mols=[1, 1], + sequence="ABA", + ) + + assert copolymer.molecules[0].n_particles == 42 + assert copolymer.molecules[1].n_particles == 62 + + def test_copolymer_with_sequence_different_num_mol( + self, polyethylene, polyDME + ): + copolymer = CoPolymer( + monomer_A=polyDME, + monomer_B=polyethylene, + lengths=[2, 3], + num_mols=[1, 2], + sequence="ABA", + ) + + assert copolymer.molecules[0].n_particles == 42 + assert copolymer.molecules[1].n_particles == 62 + assert copolymer.molecules[2].n_particles == 62 + + def test_copolymer_random_sequence(self, polyethylene, polyDME): + copolymer = CoPolymer( + monomer_A=polyDME, + monomer_B=polyethylene, + lengths=[3], + num_mols=[1], + random_sequence=True, + seed=42, + ) + # sequence is BAA + assert copolymer.n_particles == 22 + assert ("O", "C", "C", "O") in copolymer.topology_information[ + "dihedral_types" + ] diff --git a/hoomd_polymers/tests/base/test_simulation.py b/hoomd_polymers/tests/base/test_simulation.py new file mode 100644 index 00000000..49762a49 --- /dev/null +++ b/hoomd_polymers/tests/base/test_simulation.py @@ -0,0 +1,238 @@ +import os +import pickle + +import hoomd +import numpy as np + +from hoomd_polymers import Simulation +from hoomd_polymers.tests import BaseTest + + +class TestSimulate(BaseTest): + def test_initialize_from_snap(self, benzene_system): + Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + + def test_no_reference_values(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + assert np.array_equal(sim.box_lengths_reduced, sim.box_lengths) + assert sim.density_reduced == sim.density + assert sim.volume_reduced == sim.volume + assert sim.mass_reduced == sim.mass + + def test_reference_values(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.reference_values = benzene_system.reference_values + assert np.isclose(float(sim.mass.value), benzene_system.mass, atol=1e-4) + assert np.allclose(benzene_system.box.lengths, sim.box_lengths.value) + + def test_NVT(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) + assert isinstance(sim.method, hoomd.md.methods.NVT) + + def test_NPT(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.run_NPT( + kT=1.0, + n_steps=500, + pressure=0.0001, + tau_kt=0.001, + tau_pressure=0.01, + ) + assert isinstance(sim.method, hoomd.md.methods.NPT) + + def test_langevin(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.run_langevin(n_steps=500, kT=1.0, alpha=0.5) + assert isinstance(sim.method, hoomd.md.methods.Langevin) + + def test_NVE(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.run_NVE(n_steps=500) + assert isinstance(sim.method, hoomd.md.methods.NVE) + + def test_displacement_cap(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.run_displacement_cap(n_steps=500, maximum_displacement=1e-4) + assert isinstance(sim.method, hoomd.md.methods.DisplacementCapped) + + def test_update_volume(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.run_update_volume( + kT=1.0, + tau_kt=0.01, + n_steps=500, + period=1, + final_box_lengths=sim.box_lengths * 0.5, + ) + + def test_update_volume_walls(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.add_walls(wall_axis=(1, 0, 0), sigma=1.0, epsilon=1.0, r_cut=1.12) + sim.run_update_volume( + kT=1.0, + tau_kt=0.01, + n_steps=500, + period=5, + final_box_lengths=sim.box_lengths * 0.5, + ) + + def test_change_methods(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0) + assert isinstance(sim.method, hoomd.md.methods.NVT) + sim.run_NPT( + kT=1.0, tau_kt=0.01, tau_pressure=0.1, pressure=0.001, n_steps=0 + ) + assert isinstance(sim.method, hoomd.md.methods.NPT) + + def test_change_dt(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0) + sim.dt = 0.003 + sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0) + assert sim.dt == 0.003 + + def test_scale_epsilon(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + epsilons = [] + for param in sim._lj_force().params: + epsilons.append(sim._lj_force().params[param]["epsilon"]) + sim.adjust_epsilon(scale_by=0.5) + epsilons_scaled = [] + for param in sim._lj_force().params: + epsilons_scaled.append(sim._lj_force().params[param]["epsilon"]) + for i, j in zip(epsilons, epsilons_scaled): + assert np.allclose(i * 0.5, j, atol=1e-3) + + def test_shift_epsilon(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + epsilons = [] + for param in sim._lj_force().params: + epsilons.append(sim._lj_force().params[param]["epsilon"]) + sim.adjust_epsilon(shift_by=1.0) + epsilons_scaled = [] + for param in sim._lj_force().params: + epsilons_scaled.append(sim._lj_force().params[param]["epsilon"]) + for i, j in zip(epsilons, epsilons_scaled): + assert np.allclose(i + 1, j, atol=1e-3) + + def test_scale_sigma(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sigmas = [] + for param in sim._lj_force().params: + sigmas.append(sim._lj_force().params[param]["sigma"]) + sim.adjust_sigma(scale_by=0.5) + sigmas_scaled = [] + for param in sim._lj_force().params: + sigmas_scaled.append(sim._lj_force().params[param]["sigma"]) + for i, j in zip(sigmas, sigmas_scaled): + assert np.allclose(i * 0.5, j, atol=1e-3) + + def test_shift_sigma(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sigmas = [] + for param in sim._lj_force().params: + sigmas.append(sim._lj_force().params[param]["sigma"]) + sim.adjust_sigma(shift_by=1.0) + sigmas_scaled = [] + for param in sim._lj_force().params: + sigmas_scaled.append(sim._lj_force().params[param]["sigma"]) + for i, j in zip(sigmas, sigmas_scaled): + assert np.allclose(i + 1, j, atol=1e-3) + + def test_remove_force(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.remove_force(sim._lj_force()) + for i in sim.forces: + assert not isinstance(i, hoomd.md.pair.LJ) + + def test_set_integrate_group(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + assert isinstance(sim.integrate_group, hoomd.filter.All) + tag_filter = hoomd.filter.Tags([0, 1, 2, 3]) + sim.integrate_group = tag_filter + assert not isinstance(sim.integrate_group, hoomd.filter.All) + sim.run_NVT(n_steps=200, kT=1.0, tau_kt=0.01) + + def test_pickle_ff(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.pickle_forcefield("forcefield.pickle") + assert os.path.isfile("forcefield.pickle") + f = open("forcefield.pickle", "rb") + hoomd_ff = pickle.load(f) + + for i, j in zip(sim.forces, hoomd_ff): + assert type(i) is type(j) + os.remove("forcefield.pickle") + + def test_save_restart_gsd(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.save_restart_gsd("restart.gsd") + assert os.path.isfile("restart.gsd") + sim.pickle_forcefield("forcefield.pickle") + f = open("forcefield.pickle", "rb") + hoomd_ff = pickle.load(f) + Simulation(initial_state="restart.gsd", forcefield=hoomd_ff) + os.remove("forcefield.pickle") + os.remove("restart.gsd") diff --git a/hoomd_polymers/tests/base/test_system.py b/hoomd_polymers/tests/base/test_system.py new file mode 100644 index 00000000..7230078e --- /dev/null +++ b/hoomd_polymers/tests/base/test_system.py @@ -0,0 +1,244 @@ +import hoomd +import numpy as np +import unyt as u + +from hoomd_polymers import Lattice, Pack +from hoomd_polymers.library import ( + GAFF, + OPLS_AA, + OPLS_AA_DIMETHYLETHER, + OPLS_AA_PPS, +) +from hoomd_polymers.tests import BaseTest + + +class TestSystem(BaseTest): + def test_single_mol_type(self, benzene_molecule): + benzene_mols = benzene_molecule(n_mols=3) + system = Pack( + molecules=[benzene_mols], + density=0.8, + r_cut=2.5, + force_field=OPLS_AA(), + auto_scale=True, + ) + assert system.n_mol_types == 1 + assert len(system.all_molecules) == len(benzene_mols.molecules) + assert len(system.hoomd_forcefield) > 0 + assert system.n_particles == system.hoomd_snapshot.particles.N + assert system.hoomd_snapshot.particles.types == ["opls_145", "opls_146"] + assert system.reference_values.keys() == {"energy", "length", "mass"} + + def test_multiple_mol_types(self, benzene_molecule, ethane_molecule): + benzene_mol = benzene_molecule(n_mols=3) + ethane_mol = ethane_molecule(n_mols=2) + system = Pack( + molecules=[benzene_mol, ethane_mol], + density=0.8, + r_cut=2.5, + force_field=OPLS_AA(), + auto_scale=True, + ) + assert system.n_mol_types == 2 + assert len(system.all_molecules) == len(benzene_mol.molecules) + len( + ethane_mol.molecules + ) + assert system.all_molecules[0].name == "0" + assert system.all_molecules[-1].name == "1" + assert len(system.hoomd_forcefield) > 0 + assert system.n_particles == system.hoomd_snapshot.particles.N + assert system.hoomd_snapshot.particles.types == [ + "opls_135", + "opls_140", + "opls_145", + "opls_146", + ] + + def test_multiple_mol_types_different_ff( + self, dimethylether_molecule, pps_molecule + ): + dimethylether_mol = dimethylether_molecule(n_mols=3) + pps_mol = pps_molecule(n_mols=2) + system = Pack( + molecules=[dimethylether_mol, pps_mol], + density=0.8, + r_cut=2.5, + force_field=[OPLS_AA_DIMETHYLETHER(), OPLS_AA_PPS()], + auto_scale=True, + ) + assert system.n_mol_types == 2 + assert len(system.all_molecules) == len( + dimethylether_mol.molecules + ) + len(pps_mol.molecules) + assert system.all_molecules[0].name == "0" + assert system.all_molecules[-1].name == "1" + assert len(system.hoomd_forcefield) > 0 + for hoomd_force in system.hoomd_forcefield: + if isinstance(hoomd_force, hoomd.md.pair.LJ): + pairs = list(hoomd_force.params.keys()) + assert ("os", "sh") in pairs + assert system.n_particles == system.hoomd_snapshot.particles.N + assert system.hoomd_snapshot.particles.types == [ + "ca", + "ct", + "ha", + "hc", + "hs", + "os", + "sh", + ] + + def test_remove_hydrogen(self, benzene_molecule): + benzene_mol = benzene_molecule(n_mols=3) + system = Pack( + molecules=[benzene_mol], + density=0.8, + r_cut=2.5, + force_field=OPLS_AA(), + auto_scale=True, + remove_hydrogens=True, + ) + assert len(system.hoomd_forcefield) > 0 + assert list(system.hoomd_forcefield[0].params.keys()) == [ + ("opls_145", "opls_145") + ] + assert ( + system.hoomd_snapshot.particles.N + == sum(mol.n_particles for mol in benzene_mol.molecules) - 3 * 6 + ) + assert system.hoomd_snapshot.particles.types == ["opls_145"] + + def test_target_box(self, benzene_molecule): + benzene_mol = benzene_molecule(n_mols=3) + low_density_system = Pack( + molecules=[benzene_mol], + density=0.1, + r_cut=2.5, + force_field=OPLS_AA(), + auto_scale=True, + ) + high_density_system = Pack( + molecules=[benzene_mol], + density=0.9, + r_cut=2.5, + force_field=OPLS_AA(), + auto_scale=True, + ) + assert all( + low_density_system.target_box > high_density_system.target_box + ) + + def test_mass(self, pps_molecule): + pps_mol = pps_molecule(n_mols=20) + system = Pack(molecules=[pps_mol], density=1.0, r_cut=2.5) + assert np.allclose( + system.mass, ((12.011 * 6) + (1.008 * 6) + 32.06) * 20, atol=1e-4 + ) + + def test_ref_length(self, polyethylene): + polyethylene = polyethylene(lengths=5, num_mols=5) + system = Pack( + molecules=[polyethylene], + force_field=[GAFF()], + density=1.0, + r_cut=2.5, + auto_scale=True, + ) + + assert np.allclose( + system.reference_length.to("angstrom").value, 3.39966951, atol=1e-3 + ) + reduced_box = system.hoomd_snapshot.configuration.box[0:3] + calc_box = reduced_box * system.reference_length.to("nm").value + assert np.allclose(calc_box[0], system.box.Lx, atol=1e-2) + assert np.allclose(calc_box[1], system.box.Ly, atol=1e-2) + assert np.allclose(calc_box[2], system.box.Lz, atol=1e-2) + + def test_ref_mass(self, polyethylene): + polyethylene = polyethylene(lengths=5, num_mols=5) + system = Pack( + molecules=[polyethylene], + force_field=[GAFF()], + density=1.0, + r_cut=2.5, + auto_scale=True, + ) + total_red_mass = sum(system.hoomd_snapshot.particles.mass) + assert np.allclose( + system.mass, + total_red_mass * system.reference_mass.to("amu").value, + atol=1e-1, + ) + + def test_ref_energy(self, polyethylene): + polyethylene = polyethylene(lengths=5, num_mols=5) + system = Pack( + molecules=[polyethylene], + force_field=[GAFF()], + density=1.0, + r_cut=2.5, + auto_scale=True, + ) + assert np.allclose( + system.reference_energy.to("kcal/mol").value, 0.1094, atol=1e-3 + ) + + def test_ref_values_no_autoscale(self, polyethylene): + polyethylene = polyethylene(lengths=5, num_mols=5) + system = Pack( + molecules=[polyethylene], + force_field=[OPLS_AA()], + density=1.0, + r_cut=2.5, + auto_scale=False, + ) + system.reference_length = 1 * u.angstrom + system.reference_energy = 1 * u.kcal / u.mol + system.reference_mass = 1 * u.amu + assert np.allclose( + system.hoomd_snapshot.configuration.box[:3], + system.gmso_system.box.lengths.to("angstrom").value, + ) + assert dict(system.hoomd_forcefield[3].params)["opls_135", "opls_135"][ + "epsilon" + ] == system.gmso_system.sites[0].atom_type.parameters["epsilon"].to( + "kcal/mol" + ) + + def test_lattice_polymer(self, polyethylene): + polyethylene = polyethylene(lengths=2, num_mols=32) + system = Lattice( + molecules=[polyethylene], + force_field=[OPLS_AA()], + density=1.0, + r_cut=2.5, + x=1, + y=1, + n=4, + auto_scale=True, + ) + + assert system.n_mol_types == 1 + assert len(system.all_molecules) == len(polyethylene.molecules) + assert len(system.hoomd_forcefield) > 0 + assert system.n_particles == system.hoomd_snapshot.particles.N + assert system.reference_values.keys() == {"energy", "length", "mass"} + # TODO: specific asserts for lattice system? + + def test_lattice_molecule(self, benzene_molecule): + benzene_mol = benzene_molecule(n_mols=32) + system = Lattice( + molecules=[benzene_mol], + force_field=OPLS_AA(), + density=1.0, + r_cut=2.5, + x=1, + y=1, + n=4, + auto_scale=True, + ) + assert system.n_mol_types == 1 + assert len(system.all_molecules) == len(benzene_mol.molecules) + assert len(system.hoomd_forcefield) > 0 + assert system.n_particles == system.hoomd_snapshot.particles.N + assert system.reference_values.keys() == {"energy", "length", "mass"} diff --git a/hoomd_polymers/tests/base_test.py b/hoomd_polymers/tests/base_test.py index bbdf2f3c..d8848a50 100644 --- a/hoomd_polymers/tests/base_test.py +++ b/hoomd_polymers/tests/base_test.py @@ -1,35 +1,232 @@ import os +import hoomd +import mbuild as mb import pytest +from gmso.external.convert_mbuild import from_mbuild -from hoomd_polymers.systems import * -from hoomd_polymers.molecules import * -from hoomd_polymers.forcefields import * +from hoomd_polymers import Molecule, Pack, Polymer +from hoomd_polymers.library import OPLS_AA + +ASSETS_DIR = os.path.join(os.path.dirname(__file__), "assets") class BaseTest: + # @pytest.fixture(autouse=True) + # def initdir(self, tmpdir): + # tmpdir.chdir() + # + # @pytest.fixture() + # def polyethylene_system(self): + # system = Pack( + # molecule=PolyEthylene, + # n_mols=5, + # mol_kwargs={"length": 5}, + # density=0.5 + # ) + # system.apply_forcefield(forcefield=GAFF(), remove_hydrogens=False) + # return system + # + + @pytest.fixture() + def benzene_smiles(self): + return "c1ccccc1" + + @pytest.fixture() + def dimethylether_smiles(self): + return "COC" + + @pytest.fixture() + def ethane_smiles(self): + return "CC" + + @pytest.fixture() + def octane_smiles(self): + return "CCCCCCCC" + + @pytest.fixture() + def pps_smiles(self): + return "c1ccc(S)cc1" + @pytest.fixture(autouse=True) - def initdir(self, tmpdir): - tmpdir.chdir() + def benzene_mb(self, benzene_smiles): + benzene = mb.load(benzene_smiles, smiles=True) + return benzene + + @pytest.fixture() + def benzene_mol2(self): + return os.path.join(ASSETS_DIR, "benzene.mol2") + + @pytest.fixture() + def benzene_gmso(self, benzene_mb): + topology = from_mbuild(benzene_mb) + topology.identify_connections() + return topology + + @pytest.fixture() + def benzene_xml(self): + return os.path.join(ASSETS_DIR, "benzene_oplsaa.xml") + + def benzene_hoomd_pair(self, include_hydrogen=True, invalid_pair=False): + cell = hoomd.md.nlist.Cell(buffer=0.4) + lj = hoomd.md.pair.LJ(nlist=cell) + if invalid_pair: + lj.params[("C", "N")] = dict(epsilon=0.35, sigma=0.29) + lj.r_cut[("C", "N")] = 2.5 + else: + lj.params[("C", "C")] = dict(epsilon=0.35, sigma=0.29) + lj.r_cut[("C", "C")] = 2.5 + if include_hydrogen: + lj.params[("C", "H")] = dict(epsilon=0.35, sigma=0.29) + lj.r_cut[("C", "H")] = 2.5 + lj.params[("H", "H")] = dict(epsilon=0.35, sigma=0.65) + lj.r_cut[("H", "H")] = 2.5 + + return lj + + def benzene_hoomd_bond(self, include_hydrogen=True): + bond = hoomd.md.bond.Harmonic() + bond.params["C-C"] = dict(k=3.0, r0=2.38) + if include_hydrogen: + bond.params["C-H"] = dict(k=3.0, r0=2.38) + return bond + + def benzene_hoomd_angle(self, include_hydrogen=True): + angle = hoomd.md.angle.Harmonic() + angle.params["C-C-C"] = dict(k=3.0, t0=0.7851) + if include_hydrogen: + angle.params["C-C-H"] = dict(k=3.0, t0=0.7851) + return angle + + def benzene_hoomd_dihedral(self, include_hydrogen=True): + harmonic = hoomd.md.dihedral.Periodic() + harmonic.params["C-C-C-C"] = dict(k=3.0, d=0, n=1) + if include_hydrogen: + harmonic.params["C-C-C-H"] = dict(k=3.0, d=0, n=1) + harmonic.params["H-C-C-H"] = dict(k=3.0, d=-1, n=3, phi0=0) + return harmonic @pytest.fixture() - def polyethylene_system(self): + def benzene_hoomd_ff(self): + def _hoomd_ff(include_hydrogen, invalid_pair=False): + pairs = self.benzene_hoomd_pair( + include_hydrogen=include_hydrogen, invalid_pair=invalid_pair + ) + bonds = self.benzene_hoomd_bond(include_hydrogen=include_hydrogen) + angles = self.benzene_hoomd_angle(include_hydrogen=include_hydrogen) + dihedrals = self.benzene_hoomd_dihedral( + include_hydrogen=include_hydrogen + ) + return [pairs, bonds, angles, dihedrals] + + return _hoomd_ff + + @pytest.fixture() + def benzene_molecule(self, benzene_smiles): + def _benzene_molecule(n_mols): + benzene = Molecule(num_mols=n_mols, smiles=benzene_smiles) + return benzene + + return _benzene_molecule + + @pytest.fixture() + def ethane_molecule(self, ethane_smiles): + def _ethane_molecule(n_mols): + ethane = Molecule(num_mols=n_mols, smiles=ethane_smiles) + return ethane + + return _ethane_molecule + + @pytest.fixture() + def pps_molecule(self, pps_smiles): + def _pps_molecule(n_mols): + pps = Molecule(num_mols=n_mols, smiles=pps_smiles) + return pps + + return _pps_molecule + + @pytest.fixture() + def dimethylether_molecule(self, dimethylether_smiles): + def _dimethylether_molecule(n_mols): + dimethylether = Molecule( + num_mols=n_mols, smiles=dimethylether_smiles + ) + return dimethylether + + return _dimethylether_molecule + + @pytest.fixture() + def polyethylene(self, ethane_smiles): + class _PolyEthylene(Polymer): + def __init__(self, lengths, num_mols, **kwargs): + smiles = ethane_smiles + bond_indices = [2, -2] + bond_length = 0.15 + bond_orientation = [None, None] + super().__init__( + lengths=lengths, + num_mols=num_mols, + smiles=smiles, + bond_indices=bond_indices, + bond_length=bond_length, + bond_orientation=bond_orientation, + **kwargs + ) + + return _PolyEthylene + + @pytest.fixture() + def polyDME(self, dimethylether_smiles): + class _PolyDME(Polymer): + def __init__(self, lengths, num_mols, **kwargs): + smiles = dimethylether_smiles + bond_indices = [3, -1] + bond_length = 0.15 + bond_orientation = [None, None] + super().__init__( + lengths=lengths, + num_mols=num_mols, + smiles=smiles, + bond_indices=bond_indices, + bond_length=bond_length, + bond_orientation=bond_orientation, + **kwargs + ) + + return _PolyDME + + @pytest.fixture() + def benzene_system(self, benzene_mb): + benzene = Molecule(num_mols=5, compound=benzene_mb) system = Pack( - molecule=PolyEthylene, - n_mols=5, - mol_kwargs={"length": 5}, - density=0.5 + molecules=[benzene], + density=0.5, + r_cut=2.5, + force_field=OPLS_AA(), + auto_scale=True, ) - system.apply_forcefield(forcefield=GAFF(), remove_hydrogens=False) return system @pytest.fixture() - def ua_polyethylene_system(self): + def polyethylene_system(self, polyethylene): + polyethylene_mol = polyethylene(num_mols=5, lengths=5) system = Pack( - molecule=PolyEthylene, - n_mols=5, - mol_kwargs={"length": 5}, - density=0.5 + molecules=polyethylene_mol, + density=0.5, + r_cut=2.5, + force_field=OPLS_AA(), + auto_scale=True, + remove_hydrogens=True, ) - system.apply_forcefield(forcefield=GAFF(), remove_hydrogens=True) return system + + # @pytest.fixture() + # def ua_polyethylene_system(self): + # system = Pack( + # molecule=PolyEthylene, + # n_mols=5, + # mol_kwargs={"length": 5}, + # density=0.5 + # ) + # system.apply_forcefield(forcefield=GAFF(), remove_hydrogens=True) + # return system diff --git a/hoomd_polymers/tests/library/test_forcefields.py b/hoomd_polymers/tests/library/test_forcefields.py new file mode 100644 index 00000000..1accfd5b --- /dev/null +++ b/hoomd_polymers/tests/library/test_forcefields.py @@ -0,0 +1,91 @@ +import os + +import hoomd + +from hoomd_polymers.library import ( + GAFF, + OPLS_AA, + OPLS_AA_BENZENE, + OPLS_AA_DIMETHYLETHER, + OPLS_AA_PPS, + BeadSpring, + FF_from_file, +) +from hoomd_polymers.tests.base_test import ASSETS_DIR + + +class TestForceFields: + def test_GAFF(self): + ff = GAFF() + assert ff.gmso_ff is not None + + def test_OPLS_AA(self): + ff = OPLS_AA() + assert ff.gmso_ff is not None + + def test_OPLS_AA_PPS(self): + ff = OPLS_AA_PPS() + assert ff.gmso_ff is not None + + def test_OPPLS_AA_BENZENE(self): + ff = OPLS_AA_BENZENE() + assert ff.gmso_ff is not None + + def test_OPPLS_AA_DIMETHYLETHER(self): + ff = OPLS_AA_DIMETHYLETHER() + assert ff.gmso_ff is not None + + def test_FF_from_file(self): + xml_file = os.path.join(ASSETS_DIR, "test_ff.xml") + ff = FF_from_file(xml_file) + assert ff.gmso_ff is not None + + def test_BeadSpring(self): + ff = BeadSpring( + r_cut=2.5, + beads={ + "A": dict(epsilon=1.0, sigma=1.0), + "B": dict(epsilon=2.0, sigma=2.0), + }, + bonds={ + "A-A": dict(r0=1.1, k=300), + "A-B": dict(r0=1.1, k=300), + }, + angles={"A-A-A": dict(t0=2.0, k=200), "A-B-A": dict(t0=2.0, k=200)}, + dihedrals={"A-A-A-A": dict(phi0=0.0, k=100, d=-1, n=1)}, + ) + + assert isinstance(ff.hoomd_forcefield[0], hoomd.md.pair.pair.LJ) + assert isinstance(ff.hoomd_forcefield[1], hoomd.md.bond.Harmonic) + assert isinstance(ff.hoomd_forcefield[2], hoomd.md.angle.Harmonic) + assert isinstance(ff.hoomd_forcefield[3], hoomd.md.dihedral.Periodic) + + pair_types = [("A", "A"), ("A", "B"), ("B", "B")] + for param in ff.hoomd_forcefield[0].params: + assert param in pair_types + if param == ("A", "A"): + assert ff.hoomd_forcefield[0].params[param]["sigma"] == 1.0 + if param == ("B", "B"): + assert ff.hoomd_forcefield[0].params[param]["epsilon"] == 2.0 + if param == ("A", "B"): + assert ff.hoomd_forcefield[0].params[param]["epsilon"] == 1.5 + + bond_types = [("A-A"), ("A-B")] + for param in ff.hoomd_forcefield[1].params: + assert param in bond_types + assert ff.hoomd_forcefield[1].params[param]["r0"] == 1.1 + assert ff.hoomd_forcefield[1].params[param]["k"] == 300 + + angle_types = [("A-A-A"), ("A-B-A")] + for param in ff.hoomd_forcefield[2].params: + assert param in angle_types + assert ff.hoomd_forcefield[2].params[param]["t0"] == 2.0 + assert ff.hoomd_forcefield[2].params[param]["k"] == 200 + + dihedral_types = [("A-A-A-A")] + for param in ff.hoomd_forcefield[3].params: + assert param in dihedral_types + assert ff.hoomd_forcefield[3].params[param]["phi0"] == 0.0 + assert ff.hoomd_forcefield[3].params[param]["k"] == 100 + assert ff.hoomd_forcefield[3].params[param]["d"] == -1 + assert ff.hoomd_forcefield[3].params[param]["n"] == 1 diff --git a/hoomd_polymers/tests/library/test_polymers.py b/hoomd_polymers/tests/library/test_polymers.py new file mode 100644 index 00000000..6049a996 --- /dev/null +++ b/hoomd_polymers/tests/library/test_polymers.py @@ -0,0 +1,101 @@ +import mbuild as mb +import pytest + +from hoomd_polymers.library import ( + PEEK, + PPS, + LJChain, + PEKK_meta, + PEKK_para, + PolyEthylene, +) + + +class TestPolymers: + def test_pps_n_particles(self): + chain = PPS(lengths=5, num_mols=1) + monomer = mb.load(chain.smiles, smiles=True) + assert chain.n_particles == (monomer.n_particles * 5) - 8 + + def test_polyethylene(self): + chain = PolyEthylene(lengths=5, num_mols=1) + monomer = mb.load(chain.smiles, smiles=True) + assert chain.n_particles == (monomer.n_particles * 5) - 8 + + def test_pekk_meta(self): + chain = PEKK_meta(lengths=5, num_mols=1) + monomer = mb.load(chain.smiles, smiles=True) + assert chain.n_particles == (monomer.n_particles * 5) - 8 + + def test_pekk_para(self): + chain = PEKK_para(lengths=5, num_mols=1) + monomer = mb.load(chain.smiles, smiles=True) + assert chain.n_particles == (monomer.n_particles * 5) - 8 + + @pytest.mark.skip() + def test_peek(self): + chain = PEEK(lengths=5, num_mols=1) + monomer = mb.load(chain.smiles, smiles=True) + assert chain.n_particles == (monomer.n_particles * 5) - 8 + + def test_lj_chain(self): + cg_chain = LJChain( + lengths=3, + num_mols=1, + bead_sequence=["A"], + bead_mass={"A": 100}, + bond_lengths={"A-A": 1.5}, + ) + assert cg_chain.n_particles == 3 + assert cg_chain.molecules[0].mass == 300 + + def test_lj_chain_sequence(self): + cg_chain = LJChain( + lengths=3, + num_mols=1, + bead_sequence=["A", "B"], + bead_mass={"A": 100, "B": 150}, + bond_lengths={"A-A": 1.5, "A-B": 1.0}, + ) + assert cg_chain.n_particles == 6 + assert cg_chain.molecules[0].mass == 300 + 450 + + def test_lj_chain_sequence_bonds(self): + LJChain( + lengths=3, + num_mols=1, + bead_sequence=["A", "B"], + bead_mass={"A": 100, "B": 150}, + bond_lengths={"A-A": 1.5, "A-B": 1.0}, + ) + + LJChain( + lengths=3, + num_mols=1, + bead_sequence=["A", "B"], + bead_mass={"A": 100, "B": 150}, + bond_lengths={"A-A": 1.5, "B-A": 1.0}, + ) + + def test_lj_chain_sequence_bad_bonds(self): + with pytest.raises(ValueError): + LJChain( + lengths=3, + num_mols=1, + bead_sequence=["A", "B"], + bead_mass={"A": 100, "B": 150}, + bond_lengths={"A-A": 1.5}, + ) + + def test_lj_chain_sequence_bad_mass(self): + with pytest.raises(ValueError): + LJChain( + lengths=3, + num_mols=1, + bead_sequence=["A", "B"], + bead_mass={"A": 100}, + bond_lengths={"A-A": 1.5}, + ) + + def test_copolymer(self): + pass diff --git a/hoomd_polymers/tests/library/test_tensile.py b/hoomd_polymers/tests/library/test_tensile.py new file mode 100644 index 00000000..a7c4d2fe --- /dev/null +++ b/hoomd_polymers/tests/library/test_tensile.py @@ -0,0 +1,6 @@ +from hoomd_polymers.tests import BaseTest + + +class TestTensileSimulation(BaseTest): + def test_tensile(self, polyethylene_system): + pass diff --git a/hoomd_polymers/tests/modules/test_weld.py b/hoomd_polymers/tests/modules/test_weld.py new file mode 100644 index 00000000..3710e734 --- /dev/null +++ b/hoomd_polymers/tests/modules/test_weld.py @@ -0,0 +1,95 @@ +import os + +import gsd.hoomd + +from hoomd_polymers import Simulation +from hoomd_polymers.modules.welding import Interface, SlabSimulation +from hoomd_polymers.tests.base_test import BaseTest + + +class TestWelding(BaseTest): + def test_interface(self, polyethylene_system): + sim = Simulation( + initial_state=polyethylene_system.hoomd_snapshot, + forcefield=polyethylene_system.hoomd_forcefield, + log_write_freq=2000 + ) + sim.add_walls(wall_axis=(1, 0, 0), sigma=1, epsilon=1, r_cut=2) + sim.run_update_volume( + n_steps=1000, + period=10, + kT=2.0, + tau_kt=0.01, + final_box_lengths=sim.box_lengths / 2, + ) + sim.save_restart_gsd() + interface = Interface( + gsd_file="restart.gsd", interface_axis="x", gap=0.1 + ) + interface_snap = interface.hoomd_snapshot + with gsd.hoomd.open("restart.gsd", "rb") as traj: + slab_snap = traj[0] + + assert interface_snap.particles.N == slab_snap.particles.N * 2 + assert interface_snap.bonds.N == slab_snap.bonds.N * 2 + assert interface_snap.bonds.M == slab_snap.bonds.M + assert interface_snap.angles.N == slab_snap.angles.N * 2 + assert interface_snap.angles.M == slab_snap.angles.M + assert interface_snap.dihedrals.N == slab_snap.dihedrals.N * 2 + assert interface_snap.dihedrals.M == slab_snap.dihedrals.M + assert interface_snap.pairs.N == slab_snap.pairs.N * 2 + assert interface_snap.pairs.M == slab_snap.pairs.M + + if os.path.isfile("restart.gsd"): + os.remove("restart.gsd") + + def test_slab_sim_xaxis(self, polyethylene_system): + sim = SlabSimulation( + initial_state=polyethylene_system.hoomd_snapshot, + forcefield=polyethylene_system.hoomd_forcefield, + log_write_freq=2000 + ) + assert sim._axis_array == (1, 0, 0) + assert sim._axis_index == 0 + sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) + + def test_slab_sim_yaxis(self, polyethylene_system): + sim = SlabSimulation( + initial_state=polyethylene_system.hoomd_snapshot, + forcefield=polyethylene_system.hoomd_forcefield, + interface_axis="y", + log_write_freq=2000 + ) + assert sim._axis_array == (0, 1, 0) + assert sim._axis_index == 1 + sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) + + def test_slab_sim_zaxis(self, polyethylene_system): + sim = SlabSimulation( + initial_state=polyethylene_system.hoomd_snapshot, + forcefield=polyethylene_system.hoomd_forcefield, + interface_axis="z", + log_write_freq=2000 + ) + assert sim._axis_array == (0, 0, 1) + assert sim._axis_index == 2 + sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) + + def test_weld_sim(self, polyethylene_system): + sim = SlabSimulation( + initial_state=polyethylene_system.hoomd_snapshot, + forcefield=polyethylene_system.hoomd_forcefield, + log_write_freq=2000 + ) + sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) + sim.save_restart_gsd() + # Create interface system from slab restart.gsd file + interface = Interface( + gsd_file="restart.gsd", interface_axis="x", gap=0.1 + ) + sim = SlabSimulation( + initial_state=interface.hoomd_snapshot, + forcefield=polyethylene_system.hoomd_forcefield, + ) + if os.path.isfile("restart.gsd"): + os.remove("restart.gsd") diff --git a/hoomd_polymers/tests/test_molecules.py b/hoomd_polymers/tests/test_molecules.py deleted file mode 100644 index e59f8475..00000000 --- a/hoomd_polymers/tests/test_molecules.py +++ /dev/null @@ -1,93 +0,0 @@ -import os -import pytest -import random - -import numpy as np -import gsd.hoomd -from hoomd_polymers.molecules import * -#from polybinder.library import ASSETS_DIR -from base_test import BaseTest - - -class TestMolecules(BaseTest): - def test_pps(self): - chain = PPS(length=5) - monomer = mb.load(chain.smiles_str, smiles=True) - assert chain.n_particles == (monomer.n_particles*5)-8 - - def test_polyethylene(self): - chain = PolyEthylene(length=5) - monomer = mb.load(chain.smiles_str, smiles=True) - assert chain.n_particles == (monomer.n_particles*5)-8 - - def test_pekk_meta(self): - chain = PEKK_meta(length=5) - monomer = mb.load(chain.smiles_str, smiles=True) - assert chain.n_particles == (monomer.n_particles*5)-8 - - def test_pekk_para(self): - chain = PEKK_para(length=5) - monomer = mb.load(chain.smiles_str, smiles=True) - assert chain.n_particles == (monomer.n_particles*5)-8 - - @pytest.mark.skip() - def test_peek(self): - chain = PEEK(length=5) - monomer = mb.load(chain.smiles_str, smiles=True) - assert chain.n_particles == (monomer.n_particles*5)-8 - - def test_lj_chain(self): - cg_chain = LJ_chain( - length=3, - bead_sequence=["A"], - bead_mass={"A": 100}, - bond_lengths={"A-A": 1.5} - ) - assert cg_chain.n_particles == 3 - assert cg_chain.mass == 300 - - def test_lj_chain_sequence(self): - cg_chain = LJ_chain( - length=3, - bead_sequence=["A", "B"], - bead_mass={"A": 100, "B": 150}, - bond_lengths={"A-A": 1.5, "A-B": 1.0} - ) - assert cg_chain.n_particles == 6 - assert cg_chain.mass == 300 + 450 - - def test_lj_chain_sequence_bonds(self): - cg_chain = LJ_chain( - length=3, - bead_sequence=["A", "B"], - bead_mass={"A": 100, "B": 150}, - bond_lengths={"A-A": 1.5, "A-B": 1.0} - ) - - cg_chain_rev = LJ_chain( - length=3, - bead_sequence=["A", "B"], - bead_mass={"A": 100, "B": 150}, - bond_lengths={"A-A": 1.5, "B-A": 1.0} - ) - - def test_lj_chain_sequence_bad_bonds(self): - with pytest.raises(ValueError): - cg_chain = LJ_chain( - length=3, - bead_sequence=["A", "B"], - bead_mass={"A": 100, "B": 150}, - bond_lengths={"A-A": 1.5} - ) - - def test_lj_chain_sequence_bad_mass(self): - with pytest.raises(ValueError): - cg_chain = LJ_chain( - length=3, - bead_sequence=["A", "B"], - bead_mass={"A": 100}, - bond_lengths={"A-A": 1.5} - ) - - def test_copolymer(self): - pass diff --git a/hoomd_polymers/tests/test_simulations.py b/hoomd_polymers/tests/test_simulations.py deleted file mode 100644 index 12cf0463..00000000 --- a/hoomd_polymers/tests/test_simulations.py +++ /dev/null @@ -1,247 +0,0 @@ -from base_test import BaseTest -import os -import pickle -import pytest - -import gsd.hoomd -import hoomd -import numpy as np - -from hoomd_polymers.sim import Simulation - - -class TestSimulate(BaseTest): - def test_initialize_from_snap(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - - def test_no_reference_values(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - assert np.array_equal(sim.box_lengths_reduced, sim.box_lengths) - assert sim.density_reduced == sim.density - assert sim.volume_reduced == sim.volume - assert sim.mass_reduced == sim.mass - - def test_reference_values(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.reference_mass = 2 - sim.reference_distance = 2 - assert np.array_equal(2*sim.box_lengths_reduced, sim.box_lengths) - assert 2*sim.mass_reduced == sim.mass - assert sim.volume_reduced*8 == sim.volume - assert sim.density_reduced == sim.density * 4 - - def test_NVT(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) - assert isinstance(sim.method, hoomd.md.methods.NVT) - - def test_NPT(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.run_NPT( - kT=1.0, - n_steps=500, - pressure=0.0001, - tau_kt=0.001, - tau_pressure=0.01 - ) - assert isinstance(sim.method, hoomd.md.methods.NPT) - - def test_langevin(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.run_langevin(n_steps=500, kT=1.0, alpha=0.5) - assert isinstance(sim.method, hoomd.md.methods.Langevin) - - def test_NVE(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.run_NVE(n_steps=500) - assert isinstance(sim.method, hoomd.md.methods.NVE) - - def test_displacement_cap(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.run_displacement_cap(n_steps=500, maximum_displacement=1e-4) - assert isinstance(sim.method, hoomd.md.methods.DisplacementCapped) - - def test_update_volume(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.run_update_volume( - kT=1.0, - tau_kt=0.01, - n_steps=500, - period=1, - final_box_lengths=sim.box_lengths*0.5 - ) - - def test_update_volume_walls(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.add_walls(wall_axis=(1,0,0), sigma=1.0, epsilon=1.0, r_cut=1.12) - sim.run_update_volume( - kT=1.0, - tau_kt=0.01, - n_steps=500, - period=5, - final_box_lengths=sim.box_lengths*0.5 - ) - - def test_change_methods(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0) - assert isinstance(sim.method, hoomd.md.methods.NVT) - sim.run_NPT( - kT=1.0, - tau_kt=0.01, - tau_pressure=0.1, - pressure=0.001, - n_steps=0 - ) - assert isinstance(sim.method, hoomd.md.methods.NPT) - - def test_change_dt(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0) - sim.dt = 0.003 - sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0) - assert sim.dt == 0.003 - - def test_scale_epsilon(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - epsilons = [] - for param in sim._lj_force().params: - epsilons.append(sim._lj_force().params[param]["epsilon"]) - sim.adjust_epsilon(scale_by=0.5) - epsilons_scaled = [] - for param in sim._lj_force().params: - epsilons_scaled.append(sim._lj_force().params[param]["epsilon"]) - for i, j in zip(epsilons, epsilons_scaled): - assert np.allclose(i*0.5, j, atol=1e-3) - - def test_shift_epsilon(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - epsilons = [] - for param in sim._lj_force().params: - epsilons.append(sim._lj_force().params[param]["epsilon"]) - sim.adjust_epsilon(shift_by=1.0) - epsilons_scaled = [] - for param in sim._lj_force().params: - epsilons_scaled.append(sim._lj_force().params[param]["epsilon"]) - for i, j in zip(epsilons, epsilons_scaled): - assert np.allclose(i+1, j, atol=1e-3) - - def test_scale_sigma(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sigmas = [] - for param in sim._lj_force().params: - sigmas.append(sim._lj_force().params[param]["sigma"]) - sim.adjust_sigma(scale_by=0.5) - sigmas_scaled = [] - for param in sim._lj_force().params: - sigmas_scaled.append(sim._lj_force().params[param]["sigma"]) - for i, j in zip(sigmas, sigmas_scaled): - assert np.allclose(i*0.5, j, atol=1e-3) - - def test_shift_sigma(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sigmas = [] - for param in sim._lj_force().params: - sigmas.append(sim._lj_force().params[param]["sigma"]) - sim.adjust_sigma(shift_by=1.0) - sigmas_scaled = [] - for param in sim._lj_force().params: - sigmas_scaled.append(sim._lj_force().params[param]["sigma"]) - for i, j in zip(sigmas, sigmas_scaled): - assert np.allclose(i+1, j, atol=1e-3) - - def test_remove_force(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.remove_force(sim._lj_force()) - for i in sim.forces: - assert not isinstance(i, hoomd.md.pair.LJ) - - def test_set_integrate_group(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - assert isinstance(sim.integrate_group, hoomd.filter.All) - tag_filter = hoomd.filter.Tags([0, 1, 2, 3]) - sim.integrate_group = tag_filter - assert not isinstance(sim.integrate_group, hoomd.filter.All) - sim.run_NVT(n_steps=200, kT=1.0, tau_kt=0.01) - - def test_pickle_ff(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.pickle_forcefield("forcefield.pickle") - assert os.path.isfile("forcefield.pickle") - f = open("forcefield.pickle", "rb") - hoomd_ff = pickle.load(f) - - for i, j in zip(sim.forces, hoomd_ff): - assert type(i) == type(j) - os.remove("forcefield.pickle") - - def test_save_restart_gsd(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.save_restart_gsd("restart.gsd") - assert os.path.isfile("restart.gsd") - sim.pickle_forcefield("forcefield.pickle") - f = open("forcefield.pickle", "rb") - hoomd_ff = pickle.load(f) - new_sim = Simulation(initial_state="restart.gsd", forcefield=hoomd_ff) - os.remove("forcefield.pickle") - os.remove("restart.gsd") diff --git a/hoomd_polymers/tests/test_systems.py b/hoomd_polymers/tests/test_systems.py deleted file mode 100644 index a284d3b5..00000000 --- a/hoomd_polymers/tests/test_systems.py +++ /dev/null @@ -1,131 +0,0 @@ -import os -import pytest -import random -import hoomd - -from hoomd_polymers.systems import * -from hoomd_polymers.molecules import * -from hoomd_polymers.forcefields import * -from base_test import BaseTest - - -class TestSystems(BaseTest): - def test_pack(self): - system = Pack(molecule=PPS, n_mols=5, mol_kwargs={"length": 5}, density=1.0) - - def test_lattice(self): - system = Lattice( - molecule=PPS, - n_mols=32, - mol_kwargs={"length": 5}, - x=1, - y=1, - n=4, - density=1.0 - ) - - def test_set_target_box(self): - pass - - def test_hoomd_ff(self): - pass - - def test_hoomd_snap(self): - system = Pack(PEKK_para, n_mols=20, mol_kwargs={"length": 1}, density=1.0) - system.apply_forcefield(forcefield=GAFF()) - assert system.hoomd_snapshot.particles.N == system.system.n_particles - - def test_mass(self): - system = Pack(PPS, n_mols=20, mol_kwargs={"length": 1}, density=1.0) - assert np.allclose(system.mass, ((12.011*6) + (1.008*6) + 32.06)*20, atol=1e-4) - - def test_box(self): - pass - - def test_density(self): - pass - - def test_ref_distance(self): - system = Pack( - molecule=PolyEthylene, - n_mols=5, - mol_kwargs={"length": 5}, - density=1.0 - ) - system.apply_forcefield(forcefield=GAFF()) - assert np.allclose(system.reference_distance.value, 3.39966951, atol=1e-3) - reduced_box = system.hoomd_snapshot.configuration.box[0:3] - calc_box = reduced_box * system.reference_distance.to("nm").value - assert np.allclose(calc_box[0], system.box.Lx, atol=1e-2) - assert np.allclose(calc_box[1], system.box.Ly, atol=1e-2) - assert np.allclose(calc_box[2], system.box.Lz, atol=1e-2) - - def test_ref_mass(self): - system = Pack( - molecule=PolyEthylene, - n_mols=5, - mol_kwargs={"length": 5}, - density=1.0 - ) - system.apply_forcefield(forcefield=GAFF()) - total_red_mass = sum(system.hoomd_snapshot.particles.mass) - assert np.allclose( - system.mass, - total_red_mass * system.reference_mass.to("amu").value, - atol=1e-1 - ) - - def test_ref_energy(self): - system = Pack( - molecule=PolyEthylene, - n_mols=5, - mol_kwargs={"length": 5}, - density=1.0 - ) - system.apply_forcefield(forcefield=GAFF()) - assert np.allclose(system.reference_energy.value, 0.1094, atol=1e-3) - - - def test_apply_forcefield(self): - system = Pack(molecule=PolyEthylene, n_mols=5, mol_kwargs={"length": 5}, density=1.0) - system.apply_forcefield(forcefield=GAFF()) - assert isinstance(system.hoomd_snapshot, hoomd.snapshot.Snapshot) - assert isinstance(system.hoomd_forcefield, list) - - def test_remove_hydrogens(self): - system = Pack( - molecule=PolyEthylene, - n_mols=5, - mol_kwargs={"length": 5}, - density=1.0 - ) - system.apply_forcefield(forcefield=OPLS_AA(), remove_hydrogens=True) - assert system.hoomd_snapshot.particles.N == 5*5*2 - assert np.allclose( - system.mass, - sum([a.mass for a in system.typed_system.atoms]), atol=1e-1 - ) - - def test_remove_charges(self): - system = Pack( - molecule=PolyEthylene, - n_mols=5, - mol_kwargs={"length": 5}, - density=1.0 - ) - system.apply_forcefield(forcefield=OPLS_AA(), remove_charges=True) - assert sum(a.charge for a in system.typed_system.atoms) == 0 - assert sum(system.hoomd_snapshot.particles.charge) == 0 - - def test_make_charge_neutral(self): - system = Pack( - molecule=PPS, - n_mols=5, - mol_kwargs={"length": 5}, - density=1.0 - ) - system.apply_forcefield(forcefield=OPLS_AA_PPS(), make_charge_neutral=True) - assert np.allclose(0, sum([a.charge for a in system.typed_system.atoms]), atol=1e-5) - - def test_scale_parameters(self): - pass diff --git a/hoomd_polymers/tests/test_weld.py b/hoomd_polymers/tests/test_weld.py deleted file mode 100644 index 77dab263..00000000 --- a/hoomd_polymers/tests/test_weld.py +++ /dev/null @@ -1,76 +0,0 @@ -from base_test import BaseTest -import os -import pytest - -import gsd.hoomd -import hoomd -import numpy as np - -from hoomd_polymers.sim import Simulation -from hoomd_polymers.modules.welding import Interface, SlabSimulation, WeldSimulation - - -class TestWelding(BaseTest): - def test_interface(self, polyethylene_system): - sim = Simulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield - ) - sim.add_walls(wall_axis=(1,0,0), sigma=1, epsilon=1, r_cut=2) - sim.run_update_volume( - n_steps=1000, - period=10, - kT=2.0, - tau_kt=0.01, - final_box_lengths=sim.box_lengths/2, - ) - sim.save_restart_gsd() - interface = Interface(gsd_file="restart.gsd", interface_axis="x", gap=0.1) - interface_snap = interface.hoomd_snapshot - with gsd.hoomd.open("restart.gsd", "rb") as traj: - slab_snap = traj[0] - - assert interface_snap.particles.N == slab_snap.particles.N * 2 - assert interface_snap.bonds.N == slab_snap.bonds.N * 2 - assert interface_snap.bonds.M == slab_snap.bonds.M - assert interface_snap.angles.N == slab_snap.angles.N * 2 - assert interface_snap.angles.M == slab_snap.angles.M - assert interface_snap.dihedrals.N == slab_snap.dihedrals.N * 2 - assert interface_snap.dihedrals.M == slab_snap.dihedrals.M - assert interface_snap.pairs.N == slab_snap.pairs.N * 2 - assert interface_snap.pairs.M == slab_snap.pairs.M - - if os.path.isfile("restart.gsd"): - os.remove("restart.gsd") - - def test_slab_sim_xaxis(self, polyethylene_system): - sim = SlabSimulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield, - ) - assert sim._axis_array == (1,0,0) - assert sim._axis_index == 0 - sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) - - def test_slab_sim_yaxis(self, polyethylene_system): - sim = SlabSimulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield, - interface_axis="y" - ) - assert sim._axis_array == (0,1,0) - assert sim._axis_index == 1 - sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) - - def test_slab_sim_zaxis(self, polyethylene_system): - sim = SlabSimulation( - initial_state=polyethylene_system.hoomd_snapshot, - forcefield=polyethylene_system.hoomd_forcefield, - interface_axis="z" - ) - assert sim._axis_array == (0,0,1) - assert sim._axis_index == 2 - sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) - - def test_weld_sim(self): - pass diff --git a/hoomd_polymers/utils/__init__.py b/hoomd_polymers/utils/__init__.py new file mode 100644 index 00000000..96a0fce6 --- /dev/null +++ b/hoomd_polymers/utils/__init__.py @@ -0,0 +1,4 @@ +from .actions import * +from .base_types import FF_Types +from .ff_utils import xml_to_gmso_ff +from .utils import check_return_iterable, scale_charges diff --git a/hoomd_polymers/sim/actions.py b/hoomd_polymers/utils/actions.py similarity index 79% rename from hoomd_polymers/sim/actions.py rename to hoomd_polymers/utils/actions.py index 9c319d6b..79961091 100644 --- a/hoomd_polymers/sim/actions.py +++ b/hoomd_polymers/utils/actions.py @@ -7,13 +7,16 @@ def __init__(self, n_steps, sim): self.n_steps = int(n_steps) self.sim = sim self.starting_step = sim.timestep - + def act(self, timestep): if timestep != 0: - tps = np.round(self.sim.tps,2) + tps = np.round(self.sim.tps, 2) current_step = self.sim.timestep - self.starting_step - eta = np.round((self.n_steps - current_step)/(60*tps), 1) - print(f"Step {current_step} of {self.n_steps}; TPS: {tps}; ETA: {eta} minutes") + eta = np.round((self.n_steps - current_step) / (60 * tps), 1) + print( + f"Step {current_step} of {self.n_steps}; TPS: {tps}; ETA: " + f"{eta} minutes" + ) class PullParticles(hoomd.custom.Action): @@ -27,8 +30,8 @@ def act(self, timestep): with self._state.cpu_local_snapshot as snap: neg_filter = snap.particles.rtag[self.neg_filter.tags] pos_filter = snap.particles.rtag[self.pos_filter.tags] - snap.particles.position[neg_filter] -= (self.shift_by*self.axis) - snap.particles.position[pos_filter] += (self.shift_by*self.axis) + snap.particles.position[neg_filter] -= self.shift_by * self.axis + snap.particles.position[pos_filter] += self.shift_by * self.axis class UpdateWalls(hoomd.custom.Action): @@ -62,4 +65,4 @@ def __init__(self, sim, scale_factor): def act(self, timestep): self.sim.adjust_sigma(shift_by=self.scale_factor) - lj_forces = self.sim._lj_force() + self.sim._lj_force() diff --git a/hoomd_polymers/utils/base_types.py b/hoomd_polymers/utils/base_types.py new file mode 100644 index 00000000..4eea254e --- /dev/null +++ b/hoomd_polymers/utils/base_types.py @@ -0,0 +1,7 @@ +class FF_Types: + opls = "opls" + pps_opls = "pps_opls" + oplsaa = "oplsaa" + gaff = "gaff" + custom = "custom" + Hoomd = "Hoomd" diff --git a/hoomd_polymers/utils/exceptions.py b/hoomd_polymers/utils/exceptions.py new file mode 100644 index 00000000..8ce79f6d --- /dev/null +++ b/hoomd_polymers/utils/exceptions.py @@ -0,0 +1,63 @@ +class MissingPotentialError(Exception): + def __init__(self, connection="", potential_class=None): + self.connection = connection + self.potential_class = potential_class + msg = self._generate_msg() + super().__init__(msg) + + def _generate_msg(self): + return ( + f"Missing {self.potential_class} potential for" + f" {self.connection} {self.potential_type}" + ) + + @property + def potential_type(self): + return None + + +class MissingPairPotentialError(MissingPotentialError): + @property + def potential_type(self): + return "pair" + + +class MissingBondPotentialError(MissingPotentialError): + @property + def potential_type(self): + return "bond" + + +class MissingAnglePotentialError(MissingPotentialError): + @property + def potential_type(self): + return "angle" + + +class MissingDihedralPotentialError(MissingPotentialError): + @property + def potential_type(self): + return "dihedral" + + +class MissingCoulombPotentialError(MissingPotentialError): + def _generate_msg(self): + return ( + f"Missing Coulomb force {self.potential_class} " + f"for electrostatic interactions." + ) + + +class MoleculeLoadError(Exception): + def __init__(self, msg): + super().__init__(msg) + + +class ReferenceUnitError(Exception): + def __init__(self, msg): + super().__init__(msg) + + +class ForceFieldError(Exception): + def __init__(self, msg): + super().__init__(msg) diff --git a/hoomd_polymers/utils/ff_utils.py b/hoomd_polymers/utils/ff_utils.py new file mode 100644 index 00000000..d165f246 --- /dev/null +++ b/hoomd_polymers/utils/ff_utils.py @@ -0,0 +1,167 @@ +import os + +import forcefield_utilities as ffutils +import hoomd +from gmso.parameterization import apply + +from hoomd_polymers.assets import FF_DIR + +from .base_types import FF_Types +from .exceptions import ( + MissingAnglePotentialError, + MissingBondPotentialError, + MissingCoulombPotentialError, + MissingDihedralPotentialError, + MissingPairPotentialError, +) + + +def ff_xml_directory(): + ff_xml_directory = dict() + for dirpath, dirnames, filenames in os.walk(FF_DIR): + for file in filenames: + if file.endswith(".xml"): + ff_xml_directory[file.split(".xml")[0]] = os.path.join( + dirpath, file + ) + return ff_xml_directory + + +def find_xml_ff(ff_source): + xml_directory = ff_xml_directory() + if os.path.isfile(ff_source): + if not ff_source.endswith(".xml"): + raise ValueError("ForceField file type must be XML.") + ff_xml_path = ff_source + ff_type = FF_Types.custom + elif not xml_directory.get(ff_source.split(".xml")[0]): + raise ValueError( + "{} forcefield is not supported. Supported XML forcefields " + "are {}".format(ff_source, list(xml_directory.keys())) + ) + else: + ff_key = ff_source.split(".xml")[0] + ff_xml_path = xml_directory.get(ff_key) + ff_type = getattr(FF_Types, ff_key) + return ff_xml_path, ff_type + + +def xml_to_gmso_ff(ff_xml): + ff_xml_path, ff_type = find_xml_ff(ff_xml) + gmso_ff = ffutils.FoyerFFs().load(ff_xml_path).to_gmso_ff() + return gmso_ff + + +def apply_xml_ff(ff_xml_path, gmso_mol): + gmso_ff = ffutils.FoyerFFs().load(ff_xml_path).to_gmso_ff() + apply(top=gmso_mol, forcefields=gmso_ff, identify_connections=True) + # TODO: Warning if any parameter is missing? + return gmso_mol + + +def _include_hydrogen(connections, hydrogen_types): + return any(p in hydrogen_types for p in connections) + + +def _validate_hoomd_ff(forcefields, topology_information, ignore_hydrogen=True): + pair_forces = [] + bond_forces = [] + angle_forces = [] + dihedral_forces = [] + coulomb_forces = [] + + for force in forcefields: + if isinstance(force, hoomd.md.pair.Pair): + pair_forces.append(force) + elif isinstance(force, hoomd.md.bond.Bond): + bond_forces.append(force) + elif isinstance(force, hoomd.md.angle.Angle): + angle_forces.append(force) + elif isinstance(force, hoomd.md.dihedral.Dihedral): + dihedral_forces.append(force) + elif isinstance(force, hoomd.md.long_range.pppm.Coulomb): + coulomb_forces.append(force) + # check if all the required forcefields are present + if topology_information["pair_types"] and not pair_forces: + raise MissingPairPotentialError( + connection=topology_information["pair_types"], + potential_class=str(hoomd.md.pair.LJ), + ) + if topology_information["bond_types"] and not bond_forces: + raise MissingBondPotentialError( + connection=topology_information["bond_types"], + potential_class=str(hoomd.md.bond.Bond), + ) + if topology_information["angle_types"] and not angle_forces: + raise MissingAnglePotentialError( + connection=topology_information["angle_types"], + potential_class=str(hoomd.md.angle.Angle), + ) + if topology_information["dihedral_types"] and not dihedral_forces: + raise MissingDihedralPotentialError( + connection=topology_information["dihedral_types"], + potential_class=str(hoomd.md.dihedral.Dihedral), + ) + if any(topology_information["particle_charge"]) and not coulomb_forces: + raise MissingCoulombPotentialError( + potential_class=str(hoomd.md.long_range.pppm.Coulomb) + ) + + for f in pair_forces: + params = list(map(list, f.params.keys())) + for pair in topology_information["pair_types"]: + pair = list(pair) + if ignore_hydrogen and _include_hydrogen( + pair, topology_information["hydrogen_types"] + ): + # ignore pair interactions that include hydrogen atoms + continue + if not (pair in params or pair[::-1] in params): + raise MissingPairPotentialError( + connection=tuple(pair), potential_class=type(f) + ) + + for f in bond_forces: + params = list(f.params.keys()) + for bond in topology_information["bond_types"]: + bond_dir1 = "-".join(bond) + bond_dir2 = "-".join(bond[::-1]) + if ignore_hydrogen and _include_hydrogen( + bond, topology_information["hydrogen_types"] + ): + # ignore bonds that include hydrogen atoms + continue + if not (bond_dir1 in params or bond_dir2 in params): + raise MissingBondPotentialError( + connection=bond_dir1, potential_class=type(f) + ) + + for f in angle_forces: + params = list(f.params.keys()) + for angle in topology_information["angle_types"]: + angle_dir1 = "-".join(angle) + angle_dir2 = "-".join(angle[::-1]) + if ignore_hydrogen and _include_hydrogen( + angle, topology_information["hydrogen_types"] + ): + # ignore angles that include hydrogen atoms + continue + if not (angle_dir1 in params or angle_dir2 in params): + raise MissingAnglePotentialError( + connection=angle_dir1, potential_class=type(f) + ) + + for f in dihedral_forces: + params = list(f.params.keys()) + for dihedral in topology_information["dihedral_types"]: + dihedral_dir1 = "-".join(dihedral) + dihedral_dir2 = "-".join(dihedral[::-1]) + if ignore_hydrogen and _include_hydrogen( + dihedral, topology_information["hydrogen_types"] + ): + # ignore dihedrals that include hydrogen atoms + continue + if not (dihedral_dir1 in params or dihedral_dir2 in params): + raise MissingDihedralPotentialError( + connection=dihedral_dir1, potential_class=type(f) + ) diff --git a/hoomd_polymers/utils.py b/hoomd_polymers/utils/utils.py similarity index 68% rename from hoomd_polymers/utils.py rename to hoomd_polymers/utils/utils.py index 44c4d780..b87c3a39 100644 --- a/hoomd_polymers/utils.py +++ b/hoomd_polymers/utils/utils.py @@ -1,19 +1,20 @@ def check_return_iterable(obj): if isinstance(obj, dict): return [obj] + if isinstance(obj, str): + return [obj] try: iter(obj) return obj - except: + except: # noqa: E722 return [obj] -def scale_charges(charges, n_particles): +def scale_charges(charges): net_charge = sum(charges) abs_charge = sum([abs(charge) for charge in charges]) - adjust = abs(net_charge) / n_particles scaled_charges = [] for idx, charge in enumerate(charges): - new_charge = charge - (abs(charge) * (net_charge/abs_charge)) + new_charge = charge - (abs(charge) * (net_charge / abs_charge)) scaled_charges.append(new_charge) return scaled_charges diff --git a/setup.py b/setup.py index eab60ab8..01b6952e 100644 --- a/setup.py +++ b/setup.py @@ -7,20 +7,21 @@ # Note: To use the 'upload' functionality of this file, you must: # $ pip install twine -import io import os import sys from shutil import rmtree -from setuptools import find_packages, setup, Command +from setuptools import Command, find_packages, setup # Package meta-data. -NAME = 'hoomd_polymers' -DESCRIPTION = 'Package making it easier to build and simulate polymers in Hoomd-Blue' -URL = 'https://github.com/chrisjonesBSU/hoomd-polymers' -EMAIL = 'chrisjones4@u.boisestate.edu' -AUTHOR = 'Chris Jones' -REQUIRES_PYTHON = '>=3.9.0' +NAME = "hoomd_polymers" +DESCRIPTION = ( + "Package making it easier to build and simulate polymers in Hoomd-Blue" +) +URL = "https://github.com/chrisjonesBSU/hoomd-polymers" +EMAIL = "chrisjones4@u.boisestate.edu" +AUTHOR = "Chris Jones" +REQUIRES_PYTHON = ">=3.9.0" VERSION = None # What packages are required for this module to be executed? @@ -35,22 +36,22 @@ # Load the package's __version__.py module as a dictionary. about = {} if not VERSION: - with open(os.path.join(here, NAME, '__version__.py')) as f: + with open(os.path.join(here, NAME, "__version__.py")) as f: exec(f.read(), about) else: - about['__version__'] = VERSION + about["__version__"] = VERSION class UploadCommand(Command): """Support setup.py upload.""" - description = 'Build and publish the package.' + description = "Build and publish the package." user_options = [] @staticmethod def status(s): """Prints things in bold.""" - print('\033[1m{0}\033[0m'.format(s)) + print("\033[1m{0}\033[0m".format(s)) def initialize_options(self): pass @@ -60,20 +61,22 @@ def finalize_options(self): def run(self): try: - self.status('Removing previous builds…') - rmtree(os.path.join(here, 'dist')) + self.status("Removing previous builds…") + rmtree(os.path.join(here, "dist")) except OSError: pass - self.status('Building Source and Wheel (universal) distribution…') - os.system('{0} setup.py sdist bdist_wheel --universal'.format(sys.executable)) + self.status("Building Source and Wheel (universal) distribution…") + os.system( + "{0} setup.py sdist bdist_wheel --universal".format(sys.executable) + ) - self.status('Uploading the package to PyPi via Twine…') - os.system('twine upload dist/*') + self.status("Uploading the package to PyPi via Twine…") + os.system("twine upload dist/*") - self.status('Pushing git tags…') - os.system('git tag v{0}'.format(about['__version__'])) - os.system('git push --tags') + self.status("Pushing git tags…") + os.system("git tag v{0}".format(about["__version__"])) + os.system("git push --tags") sys.exit() @@ -81,41 +84,47 @@ def run(self): # Where the magic happens: setup( name=NAME, - version=about['__version__'], + version=about["__version__"], description=DESCRIPTION, author=AUTHOR, author_email=EMAIL, python_requires=REQUIRES_PYTHON, url=URL, - packages=find_packages(exclude=('tests', 'docs',)), + packages=find_packages( + exclude=( + "tests", + "docs", + ) + ), # If your package is a single module, use this instead of 'packages': # py_modules=['mypackage'], - # entry_points={ # 'console_scripts': ['mycli=mymodule:cli'], # }, - package_data={"hoomd_polymers": [ - "modules/*", - "sim/*", - "library/*", - "library/forcefields/*", - "library/monomers/*", - ]}, + package_data={ + "hoomd_polymers": [ + "modules/*", + "sim/*", + "library/*", + "library/forcefields/*", + "library/monomers/*", + ] + }, install_requires=REQUIRED, include_package_data=True, - license='MIT', + license="MIT", classifiers=[ # Trove classifiers # Full list: https://pypi.python.org/pypi?%3Aaction=list_classifiers - 'License :: OSI Approved :: GNU General Public License v3 (GPLv3)', - 'Programming Language :: Python', - 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.6', - 'Programming Language :: Python :: Implementation :: CPython', - 'Programming Language :: Python :: Implementation :: PyPy' + "License :: OSI Approved :: GNU General Public License v3 (GPLv3)", + "Programming Language :: Python", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.6", + "Programming Language :: Python :: Implementation :: CPython", + "Programming Language :: Python :: Implementation :: PyPy", ], # $ setup.py publish support. cmdclass={ - 'upload': UploadCommand, + "upload": UploadCommand, }, ) diff --git a/tutorials/Run_MD_Simulation.ipynb b/tutorials/Run_MD_Simulation.ipynb new file mode 100644 index 00000000..1f273487 --- /dev/null +++ b/tutorials/Run_MD_Simulation.ipynb @@ -0,0 +1,762 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Running a Molecular Dynamics Simulation with Hoomd-Polymers" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Overview:\n", + "In this tutorial, we will run a molecular dynamics simulation of Polyphenylene sulfide molecules using the Hoomd-Polymers package. We will use the [`HOOMD-blue`](https://hoomd-blue.readthedocs.io/en/v4.1.0/) simulation engine to run the simulation and the `hoomd_polymers` package to initialize the system of polymer configuration.\n", + "\n", + "In summary, the `hoomd-polymers` package has three main classes:\n", + "\n", + "- `Molecule`: This class is used to define the structure of a molecule (for example the structure of a polymer built from a monomer).\n", + "\n", + "- `System`: This class is used to define the system of molecules (for example a system of polymers) in a box and creates the initial `gsd` snapshot of the system. It also applies the forcefiled to the system and prepares the required forces for the simulation.\n", + "\n", + "- `Simulation`: This class is used to run the simulation using the `HOOMD-blue` simulation engine. In order to initialize a simulation, a `gsd` snapshot of the system and a list of `Hoomd` forces are required." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "import warnings\n", + "warnings.filterwarnings('ignore')" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "from hoomd_polymers.library import PPS, OPLS_AA_PPS\n", + "from hoomd_polymers import Pack, Simulation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Step 1: Defining the Molecule\n", + "In this example, we are using the pre-defined PPS molecule defined in the `hoomd_polymers` library. The `PPS` class is a subclass of the `Molecule` class. This class includes all the necessary information about the PPS molecule, including the monomer structure and how the monomers bond to create a chain. All we need to specify is the polymer length and how many polymer chain we want to create. In this example, we will create a system of 3 PPS chains with a length of 10 monomers.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "molecule = PPS(num_mols=3, lengths=10)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Step 2: Defining the System\n", + "In this step, we will use the `Pack` class, which is a subclass of the `System` class to pack a box of PPS molecules given a density. The system class creates the box and fill it with molecules, applies the force-field (if provided) to the system and creates\n", + "the initial state of the system in form a `gsd` snapshot. If force-field is provided, this class also gets the list of forces that defines the bonded and non-bonded interactions between the particles.\n", + "\n", + "In this example, we pass the molecule object created in step 1 to pack a box with density=0.8. For the force-field, we use the pre-defined `OPLS` force-field class which includes all the parameters found in the OPLS xml force-field file." + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "system = Pack(molecules=molecule, force_field=OPLS_AA_PPS(), density=0.8, r_cut=2.5, auto_scale=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can obtain the `gsd` snapshot of the system by calling the `system.snapshot` attribute.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "system.hoomd_snapshot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can also obtain the list of forces applied to the system by calling the `system.forces` attribute." + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ]" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "system.hoomd_forcefield" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Step 3: Running the Simulation\n", + "\n", + "Using the snapshot and forces provided by the system class, we can initialize the simulation. The `Simulation` class logs snapshots of the simulation in form of a `gsd` trajectory file while running simulation. The `gsd_write_freq` specifies the frequency of saving snapshots into the gsd file. This class also logs other simulation data such as timestep, potential energy, kinetic temperature, pressure and volume into a text file. The frequency for logging these information can be set by `log_write_freq` parameter." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Initializing simulation state from a snapshot.\n" + ] + } + ], + "source": [ + "sim = Simulation(initial_state=system.hoomd_snapshot, forcefield=system.hoomd_forcefield, gsd_write_freq=100, log_write_freq=100)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now run the simulation for 1000 time steps using the NVT ensemble at a given scaled temperature of 1.0." + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Step 0 of 1000; TPS: 0.0; ETA: inf minutes\n", + "Step 100 of 1000; TPS: 1100.57; ETA: 0.0 minutes\n", + "Step 200 of 1000; TPS: 1164.25; ETA: 0.0 minutes\n", + "Step 300 of 1000; TPS: 1189.75; ETA: 0.0 minutes\n", + "Step 400 of 1000; TPS: 1199.7; ETA: 0.0 minutes\n", + "Step 500 of 1000; TPS: 1204.66; ETA: 0.0 minutes\n", + "Step 600 of 1000; TPS: 1208.0; ETA: 0.0 minutes\n", + "Step 700 of 1000; TPS: 1212.62; ETA: 0.0 minutes\n", + "Step 800 of 1000; TPS: 1218.05; ETA: 0.0 minutes\n", + "Step 900 of 1000; TPS: 1218.94; ETA: 0.0 minutes\n" + ] + } + ], + "source": [ + "sim.run_NVT(n_steps=1000, kT=1.0, tau_kt=0.01)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The simulation class also allows user to run the simulation under different conditions such as NPT ensemble, NVE ensemble, Langevin dynamics. Checkout `hoomd_polymers/base/simulation.py` for more functionalities." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the rest of this tutorial, we will go through some of the features that are available in the `hoomd_polymers` package that can be tailored to specific needs." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Defining your own Molecule" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can define your own molecule in a couple of different ways:\n", + "- Using the SMILES string of the molecule\n", + "- Using the molecule file (accepted formats are: `.mol` and `.sdf`)\n", + "- Using a [`mbuild`](https://mbuild.mosdef.org/en/stable/) compound or a [`gmso`](https://gmso.mosdef.org/en/stable/) topology\n", + "- Define a subclass of the `Molecule` class" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Option 1: Using the SMILES string of the molecule" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# example of loading a molecule using the SMILES string\n", + "from hoomd_polymers import Molecule\n", + "\n", + "benzoic_acid_mol = Molecule(num_mols=20, smiles=\"c1cc(C(O)=O)ccc1\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will use `mbuild` visualization function to visualize one of the 20 benzoic acid molecules." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", + "text/html": [ + "
\n", + "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", + " jupyter labextension install jupyterlab_3dmol

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "benzoic_acid_mol.molecules[0].visualize()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Option 2: Using the molecule file" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [], + "source": [ + "# example of loading a molecule using the molecule file\n", + "\n", + "phenol_mol = Molecule(num_mols=20, file=\"../hoomd_polymers/assets/molecule_files/IPH.mol2\")" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", + "text/html": [ + "
\n", + "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", + " jupyter labextension install jupyterlab_3dmol

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "phenol_mol.molecules[0].visualize()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Option 3: Using a [`mbuild`](https://mbuild.mosdef.org/en/stable/) compound or a [`gmso`](https://gmso.mosdef.org/en/stable/) topology" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [], + "source": [ + "# example of loading a molecule from mbuild compound or gmso topology\n", + "import mbuild as mb\n", + "\n", + "mb_compound = mb.load(\"C1CCCCC1\", smiles=True)\n", + "\n", + "gmso_top = mb_compound.to_gmso()\n", + "\n", + "benzene_mol = Molecule(num_mols=20, compound=mb_compound)\n", + "benzene_mol = Molecule(num_mols=20, compound=gmso_top)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Option 4: Define a subclass of the `Molecule` class" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "is_executing": true + } + }, + "source": [ + "Checkout some examples of polymer classes defined in `hoomd_polymers/library/polymers.py`." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Defining your own Forcefield\n", + "`hoomd-polymers` package has a list of pre-defined force-fields that can be used to initialize the system. If you have the `xml` file of the forcefield, you can use the `FF_from_file` class from `hoomd_polymers.library` to create a force-field object.\n", + "You can also define your own forcefield by creating a subclass of the `foyer.Forcefield` class.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [], + "source": [ + "# example of defining a force-field using the xml file\n", + "from hoomd_polymers.library import FF_from_file\n", + "\n", + "benzene_ff = FF_from_file(xml_file=\"../hoomd_polymers/assets/forcefields/benzene_opls.xml\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Checkout `hoomd_polymers/library/forcefields.py` for more some examples of defining a forcefield using a subclass of `foyer.Forcefield` for specific molecules." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Defining your own System\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`hoomd_polymers` package has two methods of filling the box built in the `System` class: `Pack` and `Lattice`. (more info about pack and lattice?). Note that the base `System` class is considered an abstract class and cannot be called directly. " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "<__array_function__ internals>:200: DeprecationWarning: Calling nonzero on 0d arrays is deprecated, as it behaves surprisingly. Use `atleast_1d(cond).nonzero()` if the old behavior was intended. If the context of this warning is of the form `arr[nonzero(cond)]`, just use `arr[cond]`.\n", + "<__array_function__ internals>:200: DeprecationWarning: Calling nonzero on 0d arrays is deprecated, as it behaves surprisingly. Use `atleast_1d(cond).nonzero()` if the old behavior was intended. If the context of this warning is of the form `arr[nonzero(cond)]`, just use `arr[cond]`.\n" + ] + } + ], + "source": [ + "# example of defining a system using the Lattice method\n", + "\n", + "from hoomd_polymers import Lattice\n", + "from hoomd_polymers.library import OPLS_AA\n", + "\n", + "benzene_mol = Molecule(num_mols=32, smiles=\"C1CCCCC1\")\n", + "\n", + "lattice = Lattice(\n", + " molecules=[benzene_mol],\n", + " force_field=OPLS_AA(),\n", + " density=1.0,\n", + " r_cut=2.5,\n", + " x=1,\n", + " y=1,\n", + " n=4,\n", + " auto_scale=True\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", + "text/html": [ + "
\n", + "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", + " jupyter labextension install jupyterlab_3dmol

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lattice.system.visualize()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can also define your own method of filling the box by creating a subclass of the `System` class. For example, one method of filling a box with two types of molecule is creating alternate layers of each molecule type." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example of a system with multiple molecule types" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The system class can take a list of different molecule types along with different forcefields. If all molecule types use the same forcefield, then you only need to pass the forcefield once." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "from hoomd_polymers.library import OPLS_AA_DIMETHYLETHER\n", + "dimethylether_mol = Molecule(num_mols=20, smiles=\"COC\")\n", + "pps_mol = PPS(num_mols=10, lengths=4)\n", + "multi_type_system = Pack(\n", + " molecules=[dimethylether_mol, pps_mol],\n", + " density=0.8,\n", + " r_cut=2.5,\n", + " force_field=[OPLS_AA_DIMETHYLETHER(), OPLS_AA_PPS()],\n", + " auto_scale=True,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", + "text/html": [ + "
\n", + "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", + " jupyter labextension install jupyterlab_3dmol

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "multi_type_system.system.visualize()" + ] + } + ], + "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.9.16" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}