Add hoomd2 support, with tutorials

.gitignore

@@ -0,0 +1,106 @@
+# Default files
+*.txt
+*.pyc
+*.dcd
+.DS_Store
+*.egg-info
+
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+.hypothesis/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# pyenv
+.python-version
+
+# celery beat schedule file
+celerybeat-schedule
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/

msibi/optimize.py

@@ -2,7 +2,6 @@
 
 import logging
 import os
-
 import numpy as np
 
 from msibi.potentials import tail_correction
@@ -52,6 +51,7 @@ class MSIBI(object):
 
     def __init__(self, rdf_cutoff, n_rdf_points, pot_cutoff=None, r_switch=None,
                  smooth_rdfs=False, max_frames=1e3):
+
         self.states = []
         self.pairs = []
         self.n_iterations = 10  # Can be overridden in optimize().
@@ -62,13 +62,14 @@ def __init__(self, rdf_cutoff, n_rdf_points, pot_cutoff=None, r_switch=None,
         self.dr = rdf_cutoff / (n_rdf_points - 1)
         self.smooth_rdfs = smooth_rdfs
         self.rdf_r_range = np.array([0.0, self.rdf_cutoff + self.dr])
-        self.rdf_n_bins = self.n_rdf_points + 1
+        self.rdf_n_bins = self.n_rdf_points
 
-        # TODO: Description of use for pot vs rdf cutoff.
+        # Sometimes the pot_cutoff and rdf_cutoff have different ranges,
+        # e.g. to look at long-range correlations
         if not pot_cutoff:
             pot_cutoff = rdf_cutoff
         self.pot_cutoff = pot_cutoff
-        # TODO: Describe why potential needs to be messed with to match the RDF.
+
         self.pot_r = np.arange(0.0, self.pot_cutoff + self.dr, self.dr)
 
         if not r_switch:
@@ -77,14 +78,54 @@ def __init__(self, rdf_cutoff, n_rdf_points, pot_cutoff=None, r_switch=None,
 
     def optimize(self, states, pairs, n_iterations=10, engine='hoomd',
                  start_iteration=0):
+        """Optimize the pair potentials
+
+        Parameters
+        ----------
+        states : array_like, len=n_states, dtype=msibi.State
+            List of states used to optimize pair potentials.
+        pairs : array_like, len=n_pairs, dtype=msibi.Pair
+            List of pairs being optimized.
+        n_iterations : int, optional
+            Number of iterations.
+        engine : str, optional
+            Engine that runs the simulations.
+        start_iteration : int, optional
+            Start optimization at start_iteration, useful for restarting.
+
+        References
+        ----------
+        Please cite the following paper:
+
+        .. [1] T.C. Moore et al., "Derivation of coarse-grained potentials via
+           multistate iterative Boltzmann inversion," Journal of Chemical
+           Physics, vol. 140, pp. 224104, 2014.
+
         """
-        """
+
+        if engine == 'hoomd':
+            try:
+                import hoomd
+                HOOMD_VERSION = 2
+            except ImportError:
+                try:
+                    import hoomd_script
+                    HOOMD_VERSION = 1
+                except ImportError:
+                    raise UnsupportedEngine
+        else:  # Cannot find hoomd
+            HOOMD_VERSION = None
+
         for pair in pairs:
             for state, data in pair.states.items():
                 if len(data['target_rdf']) != self.n_rdf_points:
                     raise ValueError('Target RDF in {} of pair {} is not the '
                                      'same length as n_rdf_points.'.format(
-                        state.name, pair.name))
+                                     state.name, pair.name))
+
+        for state in states:
+            state.HOOMD_VERSION = HOOMD_VERSION
+
         self.states = states
         self.pairs = pairs
         self.n_iterations = n_iterations
@@ -116,18 +157,23 @@ def _recompute_rdfs(self, pair, iteration):
                          pair.states[state]['f_fit'][iteration]))
 
     def initialize(self, engine='hoomd', potentials_dir=None):
-        """Create initial table potentials and the simulation input scripts.
+        """
+        Create initial table potentials and the simulation input scripts.
 
         Parameters
         ----------
         engine : str, optional, default='hoomd'
+            Engine used to run simulations
         potentials_dir : path, optional, default="'working_dir'/potentials"
+            Directory to store potential files
 
         """
+
         if not potentials_dir:
             self.potentials_dir = os.path.join(os.getcwd(), 'potentials')
         else:
             self.potentials_dir = potentials_dir
+
         if not os.path.isdir(self.potentials_dir):
             os.mkdir(self.potentials_dir)
 
@@ -143,9 +189,11 @@ def initialize(self, engine='hoomd', potentials_dir=None):
 
             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)

msibi/state.py

@@ -3,16 +3,32 @@
 import mdtraj as md
 
 
-HOOMD_HEADER = """
+HOOMD1_HEADER = """
 from hoomd_script import *
 
 system = init.read_xml(filename="{0}", wrap_coordinates=True)
+
 T_final = {1:.1f}
 
 pot_width = {2:d}
 table = pair.table(width=pot_width)
+"""
+
+HOOMD2_HEADER = """
+import hoomd
+import hoomd.md
+from hoomd.deprecated.init import read_xml
+
+hoomd.context.initialize("")
+system = read_xml(filename="{0}", wrap_coordinates=True)
+T_final = {1:.1f}
+
+pot_width = {2:d}
+nl = hoomd.md.nlist.cell()
+table = hoomd.md.pair.table(width=pot_width, nlist=nl)
 
 """
+
 HOOMD_TABLE_ENTRY = """
 table.set_from_file('{type1}', '{type2}', filename='{potential_file}')
 """
@@ -33,22 +49,23 @@ class State(object):
         True if each query trajectory is backed up (default=False)
 
     """
-    def __init__(self, k, T, state_dir='', traj_file=None, top_file=None,
+
+    def __init__(self, kT, state_dir='', traj_file=None, top_file=None,
                  name=None, backup_trajectory=False):
-        self.kT = k * T
+        self.kT = kT
         self.state_dir = state_dir
 
         if not traj_file:
             self.traj_path = os.path.join(state_dir, 'query.dcd')
         if top_file:
             self.top_path = os.path.join(state_dir, top_file)
 
-        self.traj = None  # Will be set after first iteration.
+        self.traj = None
         if not name:
             name = 'state-{0:.3f}'.format(self.kT)
         self.name = name
 
-        self.backup_trajectory = backup_trajectory  # save trajectories?
+        self.backup_trajectory = backup_trajectory
 
     def reload_query_trajectory(self):
         """Reload the query trajectory. """
@@ -60,10 +77,15 @@ def reload_query_trajectory(self):
     def save_runscript(self, table_potentials, table_width, engine='hoomd',
                        runscript='hoomd_run_template.py'):
         """Save the input script for the MD engine. """
-
-        # TODO: Factor out for separate engines.
+
         header = list()
-        header.append(HOOMD_HEADER.format('start.xml', self.kT, table_width))
+
+        if self.HOOMD_VERSION == 1:
+            HOOMD_HEADER = HOOMD1_HEADER
+        elif self.HOOMD_VERSION == 2:
+            HOOMD_HEADER = HOOMD2_HEADER
+
+        header.append(HOOMD_HEADER.format('start.hoomdxml', self.kT, table_width))
         for type1, type2, potential_file in table_potentials:
             header.append(HOOMD_TABLE_ENTRY.format(**locals()))
         header = ''.join(header)

msibi/tutorials/lj/opt.py

@@ -1,27 +1,26 @@
 import itertools
-
+import os
 
 import numpy as np
-import os
 
 from msibi import MSIBI, State, Pair, mie
 
 
-# Clear out the temp files
 os.system('rm state*/_* rdfs/pair* potentials/* f_fits.log state*/log.txt')
 os.system('rm state*/err.txt')
+os.system('rm state*/query.dcd')
 
 # Set up global parameters.
 rdf_cutoff = 5.0
 opt = MSIBI(rdf_cutoff=rdf_cutoff, n_rdf_points=101, pot_cutoff=3.0,
         smooth_rdfs=True)
 
 # Specify states.
-state0 = State(k=1, T=0.5, state_dir='./state0', top_file='target.pdb',
+state0 = State(kT=0.5, state_dir='./state0', top_file='start.hoomdxml',
                name='state0', backup_trajectory=True)
-state1 = State(k=1, T=1.5, state_dir='./state1', top_file='target.pdb',
+state1 = State(kT=1.5, state_dir='./state1', top_file='start.hoomdxml',
                name='state1', backup_trajectory=True)
-state2 = State(k=1, T=2.0, state_dir='./state2', top_file='target.pdb',
+state2 = State(kT=2.0, state_dir='./state2', top_file='start.hoomdxml',
                name='state2', backup_trajectory=True)
 states = [state0, state1, state2]
 

msibi/tutorials/lj/state0/hoomd_run_template.py

@@ -1,14 +1,7 @@
-all = group.all()
-nvt_int = integrate.bdnvt(group=all, T=T_final)
-integrate.mode_standard(dt=0.001)
+all = hoomd.group.all()
+nvt_int = hoomd.md.integrate.langevin(group=all, kT=T_final, seed=1)
+hoomd.md.integrate.mode_standard(dt=0.001)
 
-
-run(1e2)
-output_dcd = dump.dcd(filename='query.dcd', period=100, overwrite=True)
-run(1e3)
-
-output_xml = dump.xml()
-output_xml.set_params(all=True)
-output_xml.write(filename='final.xml')
-output_pdb = dump.pdb()
-output_pdb.write(filename='final.pdb')
+hoomd.run(1e2)
+output_dcd = hoomd.dump.dcd(filename='query.dcd', period=100, overwrite=True)
+hoomd.run(1e4)

msibi/tutorials/lj/state0/start.xml → msibi/tutorials/lj/state0/start.hoomdxml

@@ -1,7 +1,7 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <hoomd_xml version="1.4">
 <configuration time_step="2100000" dimensions="3" natoms="1468" >
-<box lx="12" ly="12" lz="12"/>
+<box lx="13" ly="13" lz="13"/>
 <position num="1468">
 -4.72176265717 5.8477602005 2.7523958683
 -3.67994141579 -5.17200660706 3.88371992111