diff --git a/.codecov.yml b/.codecov.yml new file mode 100644 index 00000000..b5bd3bf7 --- /dev/null +++ b/.codecov.yml @@ -0,0 +1,11 @@ +comment: + layout: "reach, diff, flags, files" + behavior: default + require_changes: false # if true: only post the comment if coverage changes + require_base: no # [yes :: must have a base report to post] + require_head: yes # [yes :: must have a head report to post] + branches: # branch names that can post comment + - "main" +ignore: + - "msibi/tests" + - "setup.py" diff --git a/.coveragerc b/.coveragerc deleted file mode 100644 index e5893b41..00000000 --- a/.coveragerc +++ /dev/null @@ -1,8 +0,0 @@ -[report] -omit = */tutorials/* -exclude_lines = - pragma: no cover - - # Don't complain if non-runnable code isn't run: - if 0: - if __name__ == .__main__.: diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml new file mode 100644 index 00000000..f60cfdab --- /dev/null +++ b/.github/workflows/CI.yaml @@ -0,0 +1,80 @@ +name: CI + +on: + push: + branches: + - "main" + pull_request: + branches: + - "main" + schedule: + - cron: "0 0 * * *" + +jobs: + test: + if: github.event.pull_request.draft == false + name: MSIBI Tests (python) + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest] + python-version: ["3.10", "3.11", "3.12"] + + defaults: + run: + shell: bash -l {0} + + steps: + - uses: actions/checkout@v4 + name: Checkout Branch / Pull Request + + - name: Install Mamba + uses: mamba-org/setup-micromamba@v1 + with: + environment-file: environment-dev.yml + create-args: >- + python=${{ matrix.python-version }} + + - name: Install Package + run: python -m pip install -e . + + - name: Test (OS -> ${{ matrix.os }} / Python -> ${{ matrix.python-version }}) + run: python -m pytest -v --cov=msibi --cov-report=xml --cov-append --cov-config=setup.cfg --color yes --pyargs msibi + + - name: Upload Coverage Report + uses: codecov/codecov-action@v4 + with: + name: MSIBI-Coverage + verbose: true + + arch-test: + if: github.event.pull_request.draft == false + name: MSIBI Tests (arch) + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [macOS-latest, macOS-13, ubuntu-latest] + python-version: ["3.12"] + + defaults: + run: + shell: bash -l {0} + + steps: + - uses: actions/checkout@v4 + name: Checkout Branch / Pull Request + + - name: Install Mamba + uses: mamba-org/setup-micromamba@v1 + with: + environment-file: environment-dev.yml + create-args: >- + python=${{ matrix.python-version }} + + - name: Install Package + run: python -m pip install -e . + + - name: Test (OS -> ${{ matrix.os }} / Python -> ${{ matrix.python-version }}) + run: python -m pytest -v --color yes --pyargs msibi diff --git a/.github/workflows/codeql.yml b/.github/workflows/codeql.yml new file mode 100644 index 00000000..eb1c94b3 --- /dev/null +++ b/.github/workflows/codeql.yml @@ -0,0 +1,41 @@ +name: "CodeQL" + +on: + push: + branches: [ "main" ] + pull_request: + branches: [ "main" ] + schedule: + - cron: "50 6 * * 0" + +jobs: + analyze: + name: Analyze + runs-on: ubuntu-latest + permissions: + actions: read + contents: read + security-events: write + + strategy: + fail-fast: false + matrix: + language: [ python ] + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Initialize CodeQL + uses: github/codeql-action/init@v2 + with: + languages: ${{ matrix.language }} + queries: +security-and-quality + + - name: Autobuild + uses: github/codeql-action/autobuild@v2 + + - name: Perform CodeQL Analysis + uses: github/codeql-action/analyze@v2 + with: + category: "/language:${{ matrix.language }}" diff --git a/.gitignore b/.gitignore index 9e322506..1a664260 100644 --- a/.gitignore +++ b/.gitignore @@ -1,9 +1,10 @@ # Default files +*.ipynb *.pyc *.dcd .DS_Store *.egg-info - +*.idea/ # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 00000000..80fc1a44 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,32 @@ +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/astral-sh/ruff-pre-commit + rev: v0.6.1 # Ruff version + hooks: + - id: ruff + args: [--fix, --extend-ignore=E203] + - id: ruff-format + args: [--line-length=80] + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.6.0 + hooks: + - id: check-yaml + - id: end-of-file-fixer + - id: trailing-whitespace + exclude: 'msibi/tests/assets/.*' + - repo: https://github.com/pycqa/isort + rev: 5.13.2 + hooks: + - id: isort + name: isort (python) + args: + [ --profile=black, --line-length=80 ] + exclude: 'msibi/tests/assets/.* ' diff --git a/.readthedocs.yml b/.readthedocs.yml new file mode 100644 index 00000000..62c8d1e4 --- /dev/null +++ b/.readthedocs.yml @@ -0,0 +1,16 @@ +version: 2 +formats: + - htmlzip + - pdf +build: + os: ubuntu-22.04 + tools: + python: "mambaforge-4.10" + +conda: + environment: docs/environment-docs.yml + +sphinx: + builder: html + configuration: docs/source/conf.py + fail_on_warning: false diff --git a/.travis.yml b/.travis.yml deleted file mode 100755 index 147ad8ba..00000000 --- a/.travis.yml +++ /dev/null @@ -1,25 +0,0 @@ -language: generic - -sudo: false - -matrix: - include: - - { os: linux, env: PYTHON_VERSION=2.7 } - - { os: linux, env: PYTHON_VERSION=3.5 } - - { os: linux, env: PYTHON_VERSION=3.6 } - - { os: osx, env: PYTHON_VERSION=2.7 } - - { os: osx, env: PYTHON_VERSION=3.5 } - -install: - - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install md5sha1sum; fi - - source devtools/travis-ci/install.sh - - conda config --set always_yes yes --set changeps1 no - - conda env create -n test-environment python=$PYTHON_VERSION -f environment.yml - - source activate test-environment - - pip install -e . - -script: - - pytest -v --cov=msibi --cov-report= --pyargs msibi - -after_success: - - coveralls diff --git a/README.md b/README.md index c2656629..b2f3b294 100644 --- a/README.md +++ b/README.md @@ -1,20 +1,114 @@ -MultiState Iterative Boltzmann Inversion + + +# MultiState Iterative Boltzmann Inversion (MS-IBI) ---------------------------------------- -[![Build Status](https://travis-ci.org/mosdef-hub/msibi.svg?branch=master)](https://travis-ci.org/mosdef-hub/msibi) -[![Coverage Status](https://coveralls.io/repos/ctk3b/msibi/badge.svg?branch=master)](https://coveralls.io/r/ctk3b/msibi?branch=master) +[![pytest](https://github.com/mosdef-hub/msibi/actions/workflows/pytest.yml/badge.svg)](https://github.com/mosdef-hub/msibi/actions/workflows/pytest.yml) +[![codecov](https://codecov.io/gh/mosdef-hub/msibi/branch/main/graph/badge.svg?token=7NFPBMBN0I)](https://codecov.io/gh/mosdef-hub/msibi) +[![Citing MSIBI](https://img.shields.io/badge/DOI-10.1063%2F1.4880555-blue.svg)](http://dx.doi.org/10.1063/1.4880555) -A package to help you manage and run pair potential optimizations using -multistate iterative Boltzmann inversion. +A package to help you manage and run coarse-grain potential optimizations using multistate iterative Boltzmann inversion. -Install from source: -```python +## Installing MSIBI + +### Install from conda-forge (Coming soon): +``` +mamba install -c conda-forge msibi +``` + +### Install from source: +```bash git clone https://github.com/mosdef-hub/msibi.git cd msibi +mamba env create -f environment.yml +mamba activate msibi pip install . ``` +## Using MSIBI +The MSIBI package is designed to be very object oriented. Any force optimization runs requires at least one `msibi.state.State` instance, `msibi.force.Force` instance and `msibi.optimize.MSIBI` instance. More state and forces can be added as needed. Multiple forces can be added with some held fixed while others are being optimized after each iteation. MSIBI is designed to allow for optimization of both intra-molecular and inter-molecular potentials. + +MSIBI uses [Hoomd-Blue](https://hoomd-blue.readthedocs.io/en/latest/) to run optimization simulations. It is not required that you be familiar with Hoomd to use MSIBI as the simulation script is automatically generated and ran. However, it is required that you pass in the choice of [Hoomd method](https://hoomd-blue.readthedocs.io/en/latest/module-md-methods.html), [Hoomd neighbor list](https://hoomd-blue.readthedocs.io/en/latest/module-md-nlist.html), and [Hoomd thermostat](https://hoomd-blue.readthedocs.io/en/latest/module-md-methods-thermostats.html) to the `msibi.optimize.MSIBI` class. Since MSIBI utilizes Hoomd-Blue, this means that MSIBI can run on GPUs, see [Hoomd's installation guide](https://hoomd-blue.readthedocs.io/en/latest/installation.html) for instructions on ensuring your environment includes a GPU build of hoomd. + +### Example: Single state, single force +- Here is an example of learning a pair potential using a single state point with only one bead type. + +```python +import hoomd +from msibi import MSIBI, State, Pair + +optimizer = MSIBI( + nlist=hoomd.md.nlist.Cell, + integrator=hoomd.md.methods.ConstantVolume, + thermostat=hoomd.md.methods.thermostats.MTTK, + dt=0.0001, + gsd_period=int(1e4) +) + +# Create a State instance, pass in a path to the target trajectory +stateA = State(name="A", kT=2.0, traj_file="stateA.gsd", alpha=1.0, n_frames=50) +optimizer.add_state(stateA) + +# Create a Pair instance to be optimized +pairAA = Pair(type1="A", type2="A", optimize=True, r_cut=3.0, nbins=100) +# Call the set_lj() method to set an initial guess potential +pairAA.set_lj(r_min=0.001, r_cut=3.0, epsilon=1.0, sigma=1.0) +optimizer.add_force(pairAA) + +# Run 20 MSIBI iterations +optimizer.run_optimization(n_steps=2e6, n_iterations=20) +pairAA.save_potential("pairAA.csv") +``` -#### Citation [![Citing MSIBI](https://img.shields.io/badge/DOI-10.1063%2F1.4880555-blue.svg)](http://dx.doi.org/10.1063/1.4880555) +### Example: Multiple states, multiple forces +- Here is an example of learning a pair potential using multiple state points and forces. +- In this example, we set fixed bond and angle potentials that are included during iteration simulations. +- The bond potential will set a fixed harmonic force, while the angle potential will be set from a table potential file. +- This illustrates a use case of stringing together multiple MSIBI optimizations. + - For example, one MSIBI optimization can be used to learn and obtain a coarse-grained angle potnetial table file which can then be set and held fixed while learning pair potentials in a subsequent MSIBI optimization. + +```python +import hoomd +from msibi import MSIBI, State, Pair, Bond, Angle + +optimizer = MSIBI( + nlist=hoomd.md.nlist.Cell, + integrator=hoomd.md.methods.ConstantVolume, + thermostat=hoomd.md.methods.thermostats.MTTK, + dt=0.0001, + gsd_period=int(1e4) +) + +# Create 3 State instances, pass in a path to the target trajectory +stateA = State(name="A", kT=2.0, traj_file="stateA.gsd", alpha=0.2, n_frames=50) +stateB = State(name="B", kT=4.0, traj_file="stateB.gsd", alpha=0.5, n_frames=50) +stateC = State(name="C", kT=6.0, traj_file="stateC.gsd", alpha=0.3, n_frames=50) +optimizer.add_state(stateA) +optimizer.add_state(stateB) +optimizer.add_state(stateC) + +# Add bond and set a harmonic force (e.g. fit to Boltzmann inverse of the distribtion) +bondAA = Bond(type1="A", type2="A", optimize=False) +bondAA.set_harmonic(r0=1.4, k=200) +optimize.add_force(bondAA) + +# Add angle and load previously obtained table potential +angleAA = Angle(type1="A", type2="A", type3="A", optimize=False) +angleAA.set_from_file("angleAA.csv") +optimize.add_force(angleAA) + +# Create a Pair instance to be optimized. +pairAA = Pair(type1="A", type2="A", optimize=True, r_cut=3.0, nbins=100) +# Call the set_lj() method to set an initial guess potential +pairAA.set_lj(r_min=0.001, r_cut=3.0, epsilon=1.0, sigma=1.0) +optimizer.add_force(pairAA) + +# Run 20 MSIBI iterations +optimizer.run_optimization(n_steps=2e6, n_iterations=20) +pairAA.save_potential("pairAA.csv") +``` + + +## Citing MSIBI Details of the underlying method and its validation can be found [here](http://dx.doi.org/10.1063/1.4880555). If you use this package, please cite the above paper. The BibTeX reference is @@ -25,10 +119,7 @@ If you use this package, please cite the above paper. The BibTeX reference is journal = "The Journal of Chemical Physics", year = "2014", volume = "140", - number = "22", - doi = "http://dx.doi.org/10.1063/1.4880555" + number = "22", + doi = "http://dx.doi.org/10.1063/1.4880555" } ``` - - - diff --git a/devtools/travis-ci/install.sh b/devtools/travis-ci/install.sh deleted file mode 100755 index 2a7c9872..00000000 --- a/devtools/travis-ci/install.sh +++ /dev/null @@ -1,17 +0,0 @@ -#!/bin/bash - -if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then MINICONDA=Miniconda3-latest-MacOSX-x86_64.sh; fi -if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then MINICONDA=Miniconda3-latest-Linux-x86_64.sh; fi - -MINICONDA_MD5=$(curl -s https://repo.continuum.io/miniconda/ | grep -A3 $MINICONDA | sed -n '4p' | sed -n 's/ *
-from msibi.optimize import MSIBI
-from msibi.pair import Pair
-from msibi.potentials import *
-from msibi.state import State
-
-__all__ = ['MSIBI', 'Pair', 'State',
-
- # Potentials.
- 'mie', 'morse']
-
-from __future__ import division
-
-import logging
-import os
-
-import matplotlib as mpl
-try: # For use on clusters where the $DISPLAY value may not be set.
- os.environ['DISPLAY']
-except KeyError:
- mpl.use('Agg')
-
-import matplotlib.pyplot as plt
-import numpy as np
-import seaborn as sns
-
-from msibi.potentials import tail_correction
-from msibi.workers import run_query_simulations
-
-
-sns.set_style('white', {'legend.frameon': True,
- 'axes.edgecolor': '0.0',
- 'axes.linewidth': 1.0,
- 'xtick.direction': 'in',
- 'ytick.direction': 'in',
- 'xtick.major.size': 4.0,
- 'ytick.major.size': 4.0})
-
-
-[docs]class MSIBI(object):
- """Management class for orchestrating an MSIBI optimization.
-
- Parameters
- ----------
- rdf_cutoff : float
- The upper cutoff value for the RDF calculation.
- n_points : int
- The number of radius values.
- pot_cutoff : float, optional, default=rdf_cutoff
- The upper cutoff value for the potential.
- r_switch : float, optional, default=pot_r[-1] - 5 * dr
- The radius after which a tail correction is applied.
-
- Attributes
- ----------
- states : list of States
- All states to be used in the optimization procedure.
- pairs : list of Pairs
- All pairs to be used in the optimization procedure.
- n_iterations : int, optional, default=10
- The number of MSIBI iterations to perform.
- rdf_cutoff : float
- The upper cutoff value for the RDF calculation.
- n_points : int
- The number of radius values.
- dr : float, default=rdf_cutoff / (n_points - 1)
- The spacing of radius values.
- pot_r : np.ndarray, shape=(int((rdf_cutoff + dr) / dr),)
- The radius values at which the potential is computed.
- pot_cutoff : float, optional, default=rdf_cutoff
- The upper cutoff value for the potential.
- r_switch : float, optional, default=pot_r[-1] - 5 * dr
- The radius after which a tail correction is applied.
-
- """
-
- def __init__(self, rdf_cutoff, n_points, pot_cutoff=None, r_switch=None,
- status_filename='f_fits.log'):
- self.states = []
- self.pairs = []
- self.n_iterations = 10
- self.rdf_cutoff = rdf_cutoff
- self.n_points = n_points
- self.dr = rdf_cutoff / (n_points - 1)
- logging.basicConfig(filename=status_filename, level=logging.INFO,
- format='%(message)s')
-
- # TODO: description of use for pot vs rdf cutoff
- if not pot_cutoff:
- pot_cutoff = rdf_cutoff
- self.pot_cutoff = pot_cutoff
- # TODO: note on why the potential needs to be messed with to match the
- # RDF
- self.pot_r = np.linspace(0.0, pot_cutoff, n_points)
-
- if not r_switch:
- r_switch = self.pot_r[-5]
- self.r_switch = r_switch
- self.logfile = open(status_filename, 'w')
-
-[docs] def optimize(self, states, pairs, n_iterations=10, engine='hoomd'):
- """
- """
- self.states = states
- self.pairs = pairs
- if n_iterations:
- self.n_iterations = n_iterations
- self.initialize(engine=engine)
- for n in range(self.n_iterations):
- run_query_simulations(self.states, engine=engine)
-
- for pair in self.pairs:
- for state in pair.states:
- r_range = np.array([0.0, self.rdf_cutoff + 2 * self.dr])
- pair.compute_current_rdf(state, r_range, self.dr)
- pair.save_current_rdf(state, iteration=n)
- logging.info('pair {0}; state {1}; iteration {2}: {3:f}'.format(
- pair.name, state.name, n,
- pair.states[state]['f_fit'][n]))
- pair.update_potential(self.pot_r, self.r_switch)
- pair.save_table_potential(self.pot_r, self.dr, iteration=n,
- engine=engine)
-
- print("Finished iteration {0}".format(n))
-
-[docs] def initialize(self, engine='hoomd', potentials_dir=None):
- """Create initial table potentials and the simulation input scripts.
-
- Parameters
- ----------
- engine : str, optional, default='hoomd'
- potentials_dir : path, optional, default="'working_dir'/potentials"
-
- """
- if not potentials_dir:
- self.potentials_dir = os.path.join(os.getcwd(), 'potentials')
- else:
- self.potentials_dir = potentials_dir
- try:
- os.mkdir(self.potentials_dir)
- except OSError:
- # TODO: warning and maybe a "make backups" feature
- pass
-
- table_potentials = []
- for pair in self.pairs:
- potential_file = os.path.join(self.potentials_dir,
- 'pot.{0}.txt'.format(pair.name))
- pair.potential_file = potential_file
-
- table_potentials.append((pair.type1, pair.type2, potential_file))
-
- V = tail_correction(self.pot_r, pair.potential, self.r_switch)
- pair.potential = V
- # This file is written for viewing of how the potential evolves.
- pair.save_table_potential(self.pot_r, self.dr, iteration=0,
- engine=engine)
- # This file is overwritten at each iteration and actually used for
- # performing the query simulations.
- pair.save_table_potential(self.pot_r, self.dr, engine=engine)
-
- for state in self.states:
- state.save_runscript(table_potentials, table_width=len(self.pot_r),
- engine=engine)
-
-[docs] def plot(self):
- """Generate plots showing the evolution of each pair potential. """
- sns.set_palette(
- sns.cubehelix_palette(self.n_iterations, start=.5, rot=-.75))
- try:
- os.mkdir('figures')
- except OSError:
- pass
- for pair in self.pairs:
- for n in range(self.n_iterations):
- filename = 'step{0:d}.{1}'.format(
- n, os.path.basename(pair.potential_file))
- potential_file = os.path.join(self.potentials_dir, filename)
- data = np.loadtxt(potential_file)
- plt.plot(data[:, 0], data[:, 1],
- linewidth=1, label='n={0:d}'.format(n))
- plt.xlabel('r')
- plt.ylabel('V(r)')
- plt.legend()
- plt.savefig('figures/{0}.pdf'.format(pair.name))
-
-import os
-
-import mdtraj as md
-import numpy as np
-
-from six import string_types
-from msibi.utils.exceptions import UnsupportedEngine
-from msibi.utils.error_calculation import calc_similarity
-from msibi.potentials import tail_correction, head_correction, alpha_array
-
-
-[docs]class Pair(object):
- """A pair interaction to be optimized.
-
- Attributes
- ----------
- name : str
- Pair name.
- pairs : array-like, shape=(n_pairs, 2), dtype=int, optional, default=None
- Each row gives the indices of two atoms representing a pair.
- potential : func
- Values of the potential at every pot_r.
-
- """
- def __init__(self, type1, type2, potential):
- self.type1 = str(type1)
- self.type2 = str(type2)
- self.name = '{0}-{1}'.format(self.type1, self.type2)
- self.potential_file = ''
- self.states = dict()
- if isinstance(potential, string_types):
- self.potential = np.loadtxt(potential)[:, 1]
- # TODO: this could be dangerous
- else:
- self.potential = potential
- self.previous_potential = None
-
-[docs] def add_state(self, state, target_rdf, alpha, pair_indices,
- alpha_form='linear'):
- """Add a state to be used in optimizing this pair.
-
- Parameters
- ----------
- state : State
- A state object.
- target_rdf : np.ndarray, shape=(n_bins, 2), dtype=float
- Coarse-grained target RDF.
- alpha : float
- The alpha value used to scale the weight of this state.
- pair_indices : array-like, shape=(n_pairs, 2), dtype=int, optional, default=None
- Each row gives the indices of two atoms representing a pair.
-
- """
- self.states[state] = {'target_rdf': target_rdf,
- 'current_rdf': None,
- 'alpha': alpha,
- 'alpha_form': alpha_form,
- 'pair_indices': pair_indices,
- 'f_fit': []}
-
-[docs] def compute_current_rdf(self, state, r_range, dr):
- """ """
- pairs = self.states[state]['pair_indices']
- # TODO: fix units
- r, g_r = md.compute_rdf(state.traj, pairs, r_range=r_range / 10,
- bin_width=dr / 10)
- r *= 10
- rdf = np.vstack((r, g_r)).T
- self.states[state]['current_rdf'] = rdf
-
- # Compute fitness function comparing the two RDFs.
- f_fit = calc_similarity(rdf[:, 1], self.states[state]['target_rdf'][:, 1])
- self.states[state]['f_fit'].append(f_fit)
-
-[docs] def save_current_rdf(self, state, iteration):
- """ """
- filename = 'rdfs/pair_{0}-state_{1}-step{2}.txt'.format(
- self.name, state.name, iteration)
- rdf = self.states[state]['current_rdf']
- np.savetxt(filename, rdf)
-
-[docs] def update_potential(self, pot_r, r_switch=None):
- """ """
- self.previous_potential = np.copy(self.potential)
- for state in self.states:
- kT = state.kT
- alpha0 = self.states[state]['alpha']
- form = self.states[state]['alpha_form']
- alpha = alpha_array(alpha0, pot_r, form=form)
-
- current_rdf = self.states[state]['current_rdf'][:, 1]
- target_rdf = self.states[state]['target_rdf'][:, 1]
-
- # For cases where rdf_cutoff != pot_cutoff, only update the
- # potential using RDF values < pot_cutoff.
- unused_rdf_vals = current_rdf.shape[0] - self.potential.shape[0]
- if unused_rdf_vals != 0:
- current_rdf = current_rdf[:-unused_rdf_vals]
- target_rdf = target_rdf[:-unused_rdf_vals]
-
- # The actual IBI step.
- self.potential += (kT * alpha * np.log(current_rdf / target_rdf) /
- len(self.states))
-
- # Apply corrections to ensure continuous, well-behaved potentials.
- V = tail_correction(pot_r, self.potential, r_switch)
- V = head_correction(pot_r, self.potential, self.previous_potential)
- self.potential = V
-
-[docs] def save_table_potential(self, r, dr, iteration=None, engine='hoomd'):
- """Save the table potential to a file usable by the MD engine.
-
- TODO: factor out for separate engines.
- """
- V = self.potential
- F = -1.0 * np.gradient(V, dr)
- data = np.vstack([r, V, F])
-
- if iteration is not None:
- assert isinstance(iteration, int)
- basename = os.path.basename(self.potential_file)
- basename = 'step{0:d}.{1}'.format(iteration, basename)
- dirname = os.path.dirname(self.potential_file)
- iteration_filename = os.path.join(dirname, basename)
-
- if engine.lower() == 'hoomd':
- np.savetxt(self.potential_file, data.T)
- if iteration is not None:
- np.savetxt(iteration_filename, data.T)
- else:
- raise UnsupportedEngine(engine)
-
-import numpy as np
-
-from msibi.utils.general import find_nearest
-
-__all__ = ['mie', 'morse']
-
-
-[docs]def mie(r, eps, sig, m=12, n=6):
- """Mie pair potential. """
- return 4 * eps * ((sig / r) ** m - (sig / r) ** n)
-
-
-[docs]def morse(r, D, alpha, r0):
- """Morse pair potential. """
- return D * (np.exp(-2 * alpha * (r - r0)) -
- 2 * np.exp(-alpha * (r - r0)))
-
-
-def tail_correction(r, V, r_switch):
- """Apply a tail correction to a potential making it go to zero smoothly.
-
- Parameters
- ----------
- r : np.ndarray, shape=(n_points,), dtype=float
- The radius values at which the potential is given.
- V : np.ndarray, shape=r.shape, dtype=float
- The potential values at each radius value.
- r_switch : float, optional, default=pot_r[-1] - 5 * dr
- The radius after which a tail correction is applied.
-
- References
- ----------
- .. [1] https://codeblue.umich.edu/hoomd-blue/doc/classhoomd__script_1_1pair_1_1pair.html
-
- """
- r_cut = r[-1]
- idx_r_switch, r_switch = find_nearest(r, r_switch)
-
- S_r = np.ones_like(r)
- r = r[idx_r_switch:]
- S_r[idx_r_switch:] = ((r_cut ** 2 - r ** 2) ** 2 *
- (r_cut ** 2 + 2 * r ** 2 - 3 * r_switch ** 2) /
- (r_cut ** 2 - r_switch ** 2) ** 3)
- return V * S_r
-
-
-def head_correction(r, V, previous_V, form='linear'):
- """Apply head correction to V making it go to a finite value at V(0).
-
- Parameters
- ----------
- r : np.ndarray, shape=(n_points,), dtype=float
- The radius values at which the potential is given.
- V : np.ndarray, shape=r.shape, dtype=float
- The potential values at each radius value.
- previous_V : np.ndarray, shape=r.shape, dtype=float
- The potential from the previous iteration.
- form : str, optional, default='linear'
- The form of the smoothing function used.
-
- """
- if form == 'linear':
- correction_function = linear_head_correction
- else:
- raise ValueError('Unsupported head correction form: "{0}"'.format(form))
-
- for i, pot_value in enumerate(V[::-1]):
- # Apply correction function because either of the following is true:
- # * both current and target RDFs are 0 --> nan values in potential.
- # * current rdf > 0, target rdf = 0 --> +inf values in potential.
- if np.isnan(pot_value) or np.isposinf(pot_value):
- last_real = V.shape[0] - i - 1
- if last_real > len(V) - 2:
- raise RuntimeError('Undefined values in tail of potential.'
- 'This probably means you need better '
- 'sampling at this state point.')
- return correction_function(r, V, last_real)
- # Retain old potential at small r because:
- # * current rdf = 0, target rdf > 0 --> -inf values in potential.
- elif np.isneginf(pot_value):
- last_neginf = V.shape[0] - i - 1
- for i, pot_value in enumerate(V[:last_neginf+1]):
- V[i] = previous_V[i]
- return V
-
-
-def linear_head_correction(r, V, cutoff):
- """Use a linear function to smoothly force V to a finite value at V(0). """
- slope = ((V[cutoff+1] - V[cutoff+2]) / (r[cutoff+1] - r[cutoff+2]))
- V[:cutoff + 1] = slope * (r[:cutoff + 1] - r[cutoff + 1]) + V[cutoff + 1]
- return V
-
-
-def alpha_array(alpha0, pot_r, form='linear'):
- """ """
- if form == 'linear':
- return alpha0 * (1.0 - pot_r / pot_r[-1])
- else:
- raise ValueError('Unsupported alpha form')
-
-import os
-
-import mdtraj as md
-
-HOOMD_HEADER = """
-from hoomd_script import *
-
-system = init.read_xml(filename="{0}")
-T_final = {1:.1f}
-
-pot_width = {2:d}
-table = pair.table(width=pot_width)
-
-"""
-HOOMD_TABLE_ENTRY = """
-table.set_from_file('{type1}', '{type2}', filename='{potential_file}')
-"""
-
-
-[docs]class State(object):
- """A single state used as part of a multistate optimization.
-
- Attributes
- ----------
- k : float
- Boltzmann's constant in specified units.
- T : float
- Temperature in kelvin.
- traj : md.Trajectory
- The trajectory associated with this state.
- backup_trajectory : bool
- True if each query trajectory is backed up (default=False)
-
- """
- def __init__(self, k, T, state_dir='', traj_file=None, top_file=None,
- name=None, backup_trajectory=False):
- self.kT = k * T
- self.state_dir = state_dir
-
- if not traj_file:
- self.traj_path = os.path.join(state_dir, 'query.dcd')
- # TODO: check if .pdb with same name exists.
- if top_file:
- self.top_path = os.path.join(state_dir, top_file)
-
- self.traj = None # Will be set after first iteration.
- if not name:
- name = 'state-{0:.3f}'.format(self.kT)
- self.name = name
-
- self.backup_trajectory = backup_trajectory # save trajectories?
-
-[docs] def reload_query_trajectory(self):
- """ """
- if self.top_path:
- self.traj = md.load(self.traj_path, top=self.top_path)
- else:
- self.traj = md.load(self.traj_path)
-
-[docs] def save_runscript(self, table_potentials, table_width, engine='hoomd',
- runscript='hoomd_run_template.py'):
- """ """
- header = list()
- header.append(HOOMD_HEADER.format('start.xml', self.kT, table_width))
- for type1, type2, potential_file in table_potentials:
- header.append(HOOMD_TABLE_ENTRY.format(**locals()))
- header = ''.join(header)
- with open(os.path.join(self.state_dir, runscript)) as fh:
- body = ''.join(fh.readlines())
-
- runscript_file = os.path.join(self.state_dir, 'run.py')
- with open(runscript_file, 'w') as fh:
- fh.write(header)
- fh.write(body)
-
-import numpy as np
-
-
-[docs]def calc_similarity(arr1, arr2):
- f_fit = np.sum(np.absolute(arr1 - arr2))
- f_fit /= np.sum((np.absolute(arr1) + np.absolute(arr2)))
- return 1.0 - f_fit
-
-SUPPORTED_ENGINES = ['hoomd']
-
-
-[docs]class UnsupportedEngine(Exception):
- def __init__(self, engine):
- message = 'Unsupported engine: "{0}". Supported engines are: {1}'.format(
- engine, ', '.join(SUPPORTED_ENGINES))
- super(UnsupportedEngine, self).__init__(message)
-
-import glob
-import os
-import shutil
-
-import numpy as np
-
-
-[docs]def find_nearest(array, target):
- """Find array component whose numeric value is closest to 'target'. """
- idx = np.abs(array - target).argmin()
- return idx, array[idx]
-
-
-def _count_backups(filename):
- """Count the number of backups of a file in a directory. """
- head, tail = os.path.split(filename)
- backup_files = ''.join(['_.*.', tail])
- return len(glob.glob(os.path.join(head, backup_files)))
-
-
-def _backup_name(filename, n_backups):
- """Return backup filename based on the number of existing backups.
-
- Parameters
- ----------
- filename : str
- Full path to file to make backup of.
- n_backups : int
- Number of existing backups.
-
- """
- head, tail = os.path.split(filename)
- new_backup = ''.join(['_.{0:d}.'.format(n_backups), tail])
- return os.path.join(head, new_backup)
-
-
-[docs]def backup_file(filename):
- """Backup a file based on the number of backups in the file's directory.
-
- Parameters
- ----------
- filename : str
- Full path to file to make backup of.
-
- """
- n_backups = _count_backups(filename)
- new_backup = _backup_name(filename, n_backups)
- shutil.copy(filename, new_backup)
-
-import multiprocessing as mp
-from multiprocessing.dummy import Pool
-import os
-from subprocess import Popen
-
-from msibi.utils.general import backup_file
-from msibi.utils.exceptions import UnsupportedEngine
-
-
-[docs]def run_query_simulations(states, engine='hoomd'):
- """Run all query simulations for a single iteration. """
- # TODO: GPU count and proper "cluster management"
- pool = Pool(mp.cpu_count())
- print("Launching {0:d} threads...".format(mp.cpu_count()))
- if engine.lower() == 'hoomd':
- worker = _hoomd_worker
- else:
- raise UnsupportedEngine(engine)
- pool.imap(worker, states)
- pool.close()
- pool.join()
-
-
-def _hoomd_worker(state):
- """Worker for managing a single HOOMD-blue simulation. """
- log_file = os.path.join(state.state_dir, 'log.txt')
- err_file = os.path.join(state.state_dir, 'err.txt')
- with open(log_file, 'w') as log, open(err_file, 'w') as err:
- proc = Popen(['hoomd', 'run.py'], cwd=state.state_dir, stdout=log,
- stderr=err, universal_newlines=True)
- print(" Launched HOOMD in {0}...".format(state.state_dir))
- proc.communicate()
- print(" Finished in {0}.".format(state.state_dir))
- _post_query(state)
-
-
-def _post_query(state):
- state.reload_query_trajectory()
- backup_file(os.path.join(state.state_dir, 'log.txt'))
- backup_file(os.path.join(state.state_dir, 'err.txt'))
- if state.backup_trajectory:
- backup_file(state.traj_path)
-
' + _('Hide Search Matches') + '
') - .appendTo($('#searchbox')); - } - }, - - /** - * init the domain index toggle buttons - */ - initIndexTable : function() { - var togglers = $('img.toggler').click(function() { - var src = $(this).attr('src'); - var idnum = $(this).attr('id').substr(7); - $('tr.cg-' + idnum).toggle(); - if (src.substr(-9) == 'minus.png') - $(this).attr('src', src.substr(0, src.length-9) + 'plus.png'); - else - $(this).attr('src', src.substr(0, src.length-8) + 'minus.png'); - }).css('display', ''); - if (DOCUMENTATION_OPTIONS.COLLAPSE_INDEX) { - togglers.click(); - } - }, - - /** - * helper function to hide the search marks again - */ - hideSearchWords : function() { - $('#searchbox .highlight-link').fadeOut(300); - $('span.highlighted').removeClass('highlighted'); - }, - - /** - * make the url absolute - */ - makeURL : function(relativeURL) { - return DOCUMENTATION_OPTIONS.URL_ROOT + '/' + relativeURL; - }, - - /** - * get the current relative url - */ - getCurrentURL : function() { - var path = document.location.pathname; - var parts = path.split(/\//); - $.each(DOCUMENTATION_OPTIONS.URL_ROOT.split(/\//), function() { - if (this == '..') - parts.pop(); - }); - var url = parts.join('/'); - return path.substring(url.lastIndexOf('/') + 1, path.length - 1); - } -}; - -// quick alias for translations -_ = Documentation.gettext; - -$(document).ready(function() { - Documentation.init(); -}); diff --git a/docs/_build/html/_static/down-pressed.png b/docs/_build/html/_static/down-pressed.png deleted file mode 100644 index 6f7ad782..00000000 Binary files a/docs/_build/html/_static/down-pressed.png and /dev/null differ diff --git a/docs/_build/html/_static/down.png b/docs/_build/html/_static/down.png deleted file mode 100644 index 3003a887..00000000 Binary files a/docs/_build/html/_static/down.png and /dev/null differ diff --git a/docs/_build/html/_static/file.png b/docs/_build/html/_static/file.png deleted file mode 100644 index d18082e3..00000000 Binary files a/docs/_build/html/_static/file.png and /dev/null differ diff --git a/docs/_build/html/_static/fonts/fontawesome-webfont.eot b/docs/_build/html/_static/fonts/fontawesome-webfont.eot deleted file mode 100644 index 7c79c6a6..00000000 Binary files a/docs/_build/html/_static/fonts/fontawesome-webfont.eot and /dev/null differ diff --git a/docs/_build/html/_static/fonts/fontawesome-webfont.svg b/docs/_build/html/_static/fonts/fontawesome-webfont.svg deleted file mode 100644 index 45fdf338..00000000 --- a/docs/_build/html/_static/fonts/fontawesome-webfont.svg +++ /dev/null @@ -1,414 +0,0 @@ - - - \ No newline at end of file diff --git a/docs/_build/html/_static/fonts/fontawesome-webfont.ttf b/docs/_build/html/_static/fonts/fontawesome-webfont.ttf deleted file mode 100644 index e89738de..00000000 Binary files a/docs/_build/html/_static/fonts/fontawesome-webfont.ttf and /dev/null differ diff --git a/docs/_build/html/_static/fonts/fontawesome-webfont.woff b/docs/_build/html/_static/fonts/fontawesome-webfont.woff deleted file mode 100644 index 8c1748aa..00000000 Binary files a/docs/_build/html/_static/fonts/fontawesome-webfont.woff and /dev/null differ diff --git a/docs/_build/html/_static/jquery.js b/docs/_build/html/_static/jquery.js deleted file mode 100644 index 83589daa..00000000 --- a/docs/_build/html/_static/jquery.js +++ /dev/null @@ -1,2 +0,0 @@ -/*! jQuery v1.8.3 jquery.com | jquery.org/license */ -(function(e,t){function _(e){var t=M[e]={};return v.each(e.split(y),function(e,n){t[n]=!0}),t}function H(e,n,r){if(r===t&&e.nodeType===1){var i="data-"+n.replace(P,"-$1").toLowerCase();r=e.getAttribute(i);if(typeof r=="string"){try{r=r==="true"?!0:r==="false"?!1:r==="null"?null:+r+""===r?+r:D.test(r)?v.parseJSON(r):r}catch(s){}v.data(e,n,r)}else r=t}return r}function B(e){var t;for(t in e){if(t==="data"&&v.isEmptyObject(e[t]))continue;if(t!=="toJSON")return!1}return!0}function et(){return!1}function tt(){return!0}function ut(e){return!e||!e.parentNode||e.parentNode.nodeType===11}function at(e,t){do e=e[t];while(e&&e.nodeType!==1);return e}function ft(e,t,n){t=t||0;if(v.isFunction(t))return v.grep(e,function(e,r){var i=!!t.call(e,r,e);return i===n});if(t.nodeType)return v.grep(e,function(e,r){return e===t===n});if(typeof t=="string"){var r=v.grep(e,function(e){return e.nodeType===1});if(it.test(t))return v.filter(t,r,!n);t=v.filter(t,r)}return v.grep(e,function(e,r){return v.inArray(e,t)>=0===n})}function lt(e){var t=ct.split("|"),n=e.createDocumentFragment();if(n.createElement)while(t.length)n.createElement(t.pop());return n}function Lt(e,t){return e.getElementsByTagName(t)[0]||e.appendChild(e.ownerDocument.createElement(t))}function At(e,t){if(t.nodeType!==1||!v.hasData(e))return;var n,r,i,s=v._data(e),o=v._data(t,s),u=s.events;if(u){delete o.handle,o.events={};for(n in u)for(r=0,i=u[n].length;r").appendTo(i.body),n=t.css("display");t.remove();if(n==="none"||n===""){Pt=i.body.appendChild(Pt||v.extend(i.createElement("iframe"),{frameBorder:0,width:0,height:0}));if(!Ht||!Pt.createElement)Ht=(Pt.contentWindow||Pt.contentDocument).document,Ht.write(""),Ht.close();t=Ht.body.appendChild(Ht.createElement(e)),n=Dt(t,"display"),i.body.removeChild(Pt)}return Wt[e]=n,n}function fn(e,t,n,r){var i;if(v.isArray(t))v.each(t,function(t,i){n||sn.test(e)?r(e,i):fn(e+"["+(typeof i=="object"?t:"")+"]",i,n,r)});else if(!n&&v.type(t)==="object")for(i in t)fn(e+"["+i+"]",t[i],n,r);else r(e,t)}function Cn(e){return function(t,n){typeof t!="string"&&(n=t,t="*");var r,i,s,o=t.toLowerCase().split(y),u=0,a=o.length;if(v.isFunction(n))for(;u)[^>]*$|#([\w\-]*)$)/,E=/^<(\w+)\s*\/?>(?:<\/\1>|)$/,S=/^[\],:{}\s]*$/,x=/(?:^|:|,)(?:\s*\[)+/g,T=/\\(?:["\\\/bfnrt]|u[\da-fA-F]{4})/g,N=/"[^"\\\r\n]*"|true|false|null|-?(?:\d\d*\.|)\d+(?:[eE][\-+]?\d+|)/g,C=/^-ms-/,k=/-([\da-z])/gi,L=function(e,t){return(t+"").toUpperCase()},A=function(){i.addEventListener?(i.removeEventListener("DOMContentLoaded",A,!1),v.ready()):i.readyState==="complete"&&(i.detachEvent("onreadystatechange",A),v.ready())},O={};v.fn=v.prototype={constructor:v,init:function(e,n,r){var s,o,u,a;if(!e)return this;if(e.nodeType)return this.context=this[0]=e,this.length=1,this;if(typeof e=="string"){e.charAt(0)==="<"&&e.charAt(e.length-1)===">"&&e.length>=3?s=[null,e,null]:s=w.exec(e);if(s&&(s[1]||!n)){if(s[1])return n=n instanceof v?n[0]:n,a=n&&n.nodeType?n.ownerDocument||n:i,e=v.parseHTML(s[1],a,!0),E.test(s[1])&&v.isPlainObject(n)&&this.attr.call(e,n,!0),v.merge(this,e);o=i.getElementById(s[2]);if(o&&o.parentNode){if(o.id!==s[2])return r.find(e);this.length=1,this[0]=o}return this.context=i,this.selector=e,this}return!n||n.jquery?(n||r).find(e):this.constructor(n).find(e)}return v.isFunction(e)?r.ready(e):(e.selector!==t&&(this.selector=e.selector,this.context=e.context),v.makeArray(e,this))},selector:"",jquery:"1.8.3",length:0,size:function(){return this.length},toArray:function(){return l.call(this)},get:function(e){return e==null?this.toArray():e<0?this[this.length+e]:this[e]},pushStack:function(e,t,n){var r=v.merge(this.constructor(),e);return r.prevObject=this,r.context=this.context,t==="find"?r.selector=this.selector+(this.selector?" ":"")+n:t&&(r.selector=this.selector+"."+t+"("+n+")"),r},each:function(e,t){return v.each(this,e,t)},ready:function(e){return v.ready.promise().done(e),this},eq:function(e){return e=+e,e===-1?this.slice(e):this.slice(e,e+1)},first:function(){return this.eq(0)},last:function(){return this.eq(-1)},slice:function(){return this.pushStack(l.apply(this,arguments),"slice",l.call(arguments).join(","))},map:function(e){return this.pushStack(v.map(this,function(t,n){return e.call(t,n,t)}))},end:function(){return this.prevObject||this.constructor(null)},push:f,sort:[].sort,splice:[].splice},v.fn.init.prototype=v.fn,v.extend=v.fn.extend=function(){var e,n,r,i,s,o,u=arguments[0]||{},a=1,f=arguments.length,l=!1;typeof u=="boolean"&&(l=u,u=arguments[1]||{},a=2),typeof u!="object"&&!v.isFunction(u)&&(u={}),f===a&&(u=this,--a);for(;at |
- m | ||
- | - msibi | - |
- | - msibi.optimize | - |
- | - msibi.pair | - |
- | - msibi.potentials | - |
- | - msibi.state | - |
- | - msibi.utils | - |
- | - msibi.utils.error_calculation | - |
- | - msibi.utils.exceptions | - |
- | - msibi.utils.general | - |
- | - msibi.workers | - |
- Please activate JavaScript to enable the search - functionality. -
-- From here you can search these documents. Enter your search - words into the box below and click "search". Note that the search - function will automatically search for all of the words. Pages - containing fewer words won't appear in the result list. -
- - -Coming soon!
-