Finish fluctmatch-3.4.0

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 3.3.5
+current_version = 3.4.0
 commit = True
 tag = True
 

CHANGELOG.rst

@@ -2,6 +2,13 @@
 Changelog
 =========
 
+3.3.5 (2018-07-09)
+------------------
+
+* splittraj now writes Gromacs log files in each data subdirectory
+* table_convert properly assigns residue names
+* Changed RTF column length to 10
+
 3.3.3 (2018-05-16)
 ------------------
 

README.rst

@@ -44,9 +44,9 @@ Overview
     :alt: PyPI Package latest release
     :target: https://pypi.python.org/pypi/fluctmatch
 
-.. |commits-since| image:: https://img.shields.io/github/commits-since/tclick/python-fluctmatch/v3.3.5.svg
+.. |commits-since| image:: https://img.shields.io/github/commits-since/tclick/python-fluctmatch/v3.4.0.svg
     :alt: Commits since latest release
-    :target: https://github.com/tclick/python-fluctmatch/compare/v3.3.5...master
+    :target: https://github.com/tclick/python-fluctmatch/compare/v3.4.0...master
 
 .. |wheel| image:: https://img.shields.io/pypi/wheel/fluctmatch.svg
     :alt: PyPI Wheel

docs/conf.py

@@ -26,7 +26,7 @@
 year = u'2014-2017'
 author = u'Timothy Click'
 copyright = '{0}, {1}'.format(year, author)
-version = release = u'3.3.5'
+version = release = u'3.4.0'
 
 pygments_style = 'trac'
 templates_path = ['.']

setup.py

@@ -45,7 +45,7 @@ def read(*names, **kwargs):
 
 setup(
     name="fluctmatch",
-    version="3.3.5",
+    version="3.4.0",
     license="BSD license",
     description="Elastic network model using fluctuation matching.",
     long_description="%s\n%s" % (

src/fluctmatch/__init__.py

@@ -24,7 +24,7 @@
 
 import logging
 
-__version__ = "3.3.5"
+__version__ = "3.4.0"
 
 _MODELS = dict()
 _DESCRIBE = dict()

src/fluctmatch/analysis/entropy.py

@@ -120,10 +120,12 @@ def windiff_entropy(self, bins=100):
         -------
         Entropy difference between two windows
         """
+
         # Calculate value of maximum probability and define penalty value.
-        header = ["segidI", "resI"]
-        normalize = lambda x: x / x.sum()
+        def normalize(x):
+            return x / x.sum()
 
+        header = ["segidI", "resI"]
         table = self._table._separate(self._table.table)
         hist, edges = np.histogram(
             table, range=(1e-4, table.values.max()), bins=bins)

src/fluctmatch/analysis/fluctsca.py

@@ -27,6 +27,7 @@
 import os
 from os import path
 
+import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
 from scipy import linalg
@@ -170,7 +171,6 @@ def figUnits(v1,
      '''
     import colorsys
 
-    Ntot = len(v1)
     # Plot all items in white:
     fig = plt.figure()
     ax = fig.add_subplot(111, projection='3d')
@@ -250,41 +250,39 @@ def icList(Vpica, kpos, Csca, p_cut=0.95):
     icList, icsize, sortedpos, cutoff, pd  = icList(Vsca,Lsca,Lrand) """
     #do the PDF/CDF fit, and assign cutoffs
     Npos = len(Vpica)
-    cutoff = []
-    scaled_pdf = []
-    all_fits = []
+    cutoff = list()
+    scaled_pdf = list()
+    all_fits = list()
     for k in range(kpos):
-        pd = t.fit(Vpica[:, k])
+        pd = t.fit(Vpica[:,k])
         all_fits.append(pd)
-        iqr = scoreatpercentile(Vpica[:, k], 75) - scoreatpercentile(
-            Vpica[:, k], 25)
-        binwidth = 2 * iqr * (len(Vpica[:, k])**(-0.33))
-        nbins = round((max(Vpica[:, k]) - min(Vpica[:, k])) / binwidth)
-        h_params = np.histogram(Vpica[:, k], nbins.astype(np.int))
+        iqr = scoreatpercentile(Vpica[:,k],75) - scoreatpercentile(Vpica[:,k],25)
+        binwidth=2*iqr*(len(Vpica[:,k])**(-0.33))
+        nbins=round((max(Vpica[:,k])-min(Vpica[:,k]))/binwidth)
+        h_params = np.histogram(Vpica[:,k], int(nbins))
         x_dist = np.linspace(min(h_params[1]), max(h_params[1]), num=100)
-        area_hist = Npos * (h_params[1][2] - h_params[1][1])
-        scaled_pdf.append(area_hist * (t.pdf(x_dist, pd[0], pd[1], pd[2])))
-        cd = t.cdf(x_dist, pd[0], pd[1], pd[2])
+        area_hist=Npos*(h_params[1][2]-h_params[1][1]);
+        scaled_pdf.append(area_hist*(t.pdf(x_dist,pd[0],pd[1],pd[2])))
+        cd = t.cdf(x_dist,pd[0],pd[1],pd[2])
         tmp = scaled_pdf[k].argmax()
-        if abs(max(Vpica[:, k])) > abs(min(Vpica[:, k])):
+        if abs(max(Vpica[:,k])) > abs(min(Vpica[:,k])):
             tail = cd[tmp:len(cd)]
         else:
             cd = 1 - cd
             tail = cd[0:tmp]
-        diff = abs(tail - p_cut)
+        diff = abs(tail - p_cut);
         x_pos = diff.argmin()
-        cutoff.append(x_dist[x_pos + tmp])
-    #select the positions with significant contributions to each IC
-    ic_init = []
-    for k in range(kpos):
-        ic_init.append([i for i in range(Npos) if Vpica[i, k] > cutoff[k]])
-    #construct the sorted, non-redundant iclist
-    sortedpos = []
-    icsize = []
-    ics = []
+        cutoff.append(x_dist[x_pos+tmp])
+    # select the positions with significant contributions to each IC
+    ic_init = list()
+    for k in range(kpos): ic_init.append([i for i in range(Npos) if Vpica[i,k]> cutoff[k]])
+    # construct the sorted, non-redundant iclist
+    sortedpos = list()
+    icsize = list()
+    ics = list()
+    icpos_tmp = list()
     Csca_nodiag = Csca.copy()
-    for i in range(Npos):
-        Csca_nodiag[i, i] = 0
+    for i in range(Npos): Csca_nodiag[i,i]=0
     for k in range(kpos):
         icpos_tmp = list(ic_init[k])
         for kprime in [kp for kp in range(kpos) if (kp != k)]:
@@ -293,12 +291,12 @@ def icList(Vpica, kpos, Csca, p_cut=0.95):
                 remsec = np.linalg.norm(Csca_nodiag[i,ic_init[k]]) \
                          < np.linalg.norm(Csca_nodiag[i,ic_init[kprime]])
                 if remsec: icpos_tmp.remove(i)
-        sortedpos += sorted(icpos_tmp, key=lambda i: -Vpica[i, k])
+        sortedpos += sorted(icpos_tmp, key=lambda i: -Vpica[i,k])
         icsize.append(len(icpos_tmp))
         s = Unit()
-        s.items = sorted(icpos_tmp, key=lambda i: -Vpica[i, k])
-        s.col = k / kpos
-        s.vect = -Vpica[s.items, k]
+        s.items = sorted(icpos_tmp, key=lambda i: -Vpica[i,k])
+        s.col = k/kpos
+        s.vect = -Vpica[s.items,k]
         ics.append(s)
     return ics, icsize, sortedpos, cutoff, scaled_pdf, all_fits
 
@@ -315,20 +313,21 @@ def basicICA(x, r, Niter):
       -  `w` = unmixing matrix
       -  `change` = record of incremental changes during the iterations.
 
-    **Note:** r and Niter should be adjusted to achieve convergence, which should be assessed by visualizing "change" with plot(range(iter) ,change)
+    **Note:** r and Niter should be adjusted to achieve convergence, which should be assessed by visualizing 'change' with plot(range(iter) ,change)
 
     **Example:**
       >>> [w, change] = basicICA(x, r, Niter)
 
     """
     [L, M] = x.shape
     w = np.eye(L)
-    change = []
+    change = list()
     for _ in range(Niter):
         w_old = np.copy(w)
         u = w.dot(x)
-        w += r * (M * np.eye(L) + (1 - 2 * (1. /
-                                            (1 + np.exp(-u)))).dot(u.T)).dot(w)
+        w += r * (
+                M * np.eye(L) + (1 - 2 * (1. / (1 + np.exp(-u)))).dot(u.T)).dot(
+            w)
         delta = (w - w_old).ravel()
         change.append(delta.dot(delta.T))
     return [w, change]

src/fluctmatch/commands/cmd_sca.py

@@ -133,8 +133,8 @@
         file_okay=True,
         resolve_path=True,
     ))
-def cli(logfile, ntrials, std, kpos, pcut, ressep, output, subset, transformation,
-        filename):
+def cli(logfile, ntrials, std, kpos, pcut, ressep, output, subset,
+        transformation, filename):
     # Setup logger
     logging.config.dictConfig({
         "version": 1,

src/fluctmatch/commands/cmd_thermo.py

@@ -159,6 +159,9 @@ def cli(datadir, logfile, outdir, topology, trajectory, nma_exec, temperature,
     except OSError:
         pass
 
+    logger.info("Calculating thermodynamic properties.")
+    logger.warning("Depending upon the size of the system, this may take a "
+                   "while.")
     thermodynamics.create_thermo_tables(
         datadir,
         outdir,

src/fluctmatch/commands/cmd_write_charmm.py

@@ -177,4 +177,5 @@ def cli(
         write_traj=write_traj,
     )
     universe = mda.Universe(topology, trajectory)
+    logger.info("Writing CHARMM files.")
     fmutils.write_charmm_files(universe, **kwargs)

src/fluctmatch/commands/cmd_write_rtf.py

@@ -0,0 +1,139 @@
+# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding: utf-8 -*-
+# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
+#
+# pysca --- https://github.com/tclick/python-pysca
+# Copyright (c) 2015-2017 The pySCA Development Team and contributors
+# (see the file AUTHORS for the full list of names)
+#
+# Released under the New BSD license.
+#
+# Please cite your use of fluctmatch in published work:
+#
+# Timothy H. Click, Nixon Raj, and Jhih-Wei Chu.
+# Calculation of Enzyme Fluctuograms from All-Atom Molecular Dynamics
+# Simulation. Meth Enzymology. 578 (2016), 327-342,
+# doi:10.1016/bs.mie.2016.05.024.
+#
+from __future__ import (
+    absolute_import,
+    division,
+    print_function,
+    unicode_literals,
+)
+
+from future.utils import (native_str)
+
+import logging
+import logging.config
+import os
+from os import path
+
+import click
+import MDAnalysis as mda
+from fluctmatch.topology import RTF
+
+
+@click.command(
+    "write_rtf", short_help="Create an RTF file from a structure file.")
+@click.option(
+    "-s",
+    "topology",
+    metavar="FILE",
+    default=path.join(os.getcwd(), "md.tpr"),
+    show_default=True,
+    type=click.Path(exists=False, file_okay=True, resolve_path=True),
+    help="Gromacs topology file (e.g., tpr gro g96 pdb brk ent)",
+)
+@click.option(
+    "-f",
+    "trajectory",
+    metavar="FILE",
+    default=path.join(os.getcwd(), "md.pdb"),
+    show_default=True,
+    type=click.Path(exists=False, file_okay=True, resolve_path=True),
+    help="Trajectory file (e.g. xtc trr dcd, crd, cor, pdb, tpr)",
+)
+@click.option(
+    "-l",
+    "--logfile",
+    metavar="LOG",
+    show_default=True,
+    default=path.join(os.getcwd(), "write_rtf.log"),
+    type=click.Path(exists=False, file_okay=True, resolve_path=True),
+    help="Log file",
+)
+@click.option(
+    "-o",
+    "--outfile",
+    metavar="FILE",
+    show_default=True,
+    default=path.join(os.getcwd(), "md.rtf"),
+    type=click.Path(exists=False, file_okay=True, resolve_path=True),
+    help="CHARMM topology file",
+)
+@click.option(
+    "--no-decl",
+    "decl",
+    is_flag=True,
+    help="Include declaration section in CHARMM topology file",
+)
+@click.option(
+    "--uniform",
+    "mass",
+    is_flag=True,
+    help="Set uniform mass of beads to 1.0",
+)
+def cli(
+        topology,
+        trajectory,
+        logfile,
+        outfile,
+        decl,
+        mass,
+):
+    logging.config.dictConfig({
+        "version": 1,
+        "disable_existing_loggers": False,  # this fixes the problem
+        "formatters": {
+            "standard": {
+                "class": "logging.Formatter",
+                "format": "%(name)-12s %(levelname)-8s %(message)s",
+            },
+            "detailed": {
+                "class": "logging.Formatter",
+                "format":
+                "%(asctime)s %(name)-15s %(levelname)-8s %(message)s",
+                "datefmt": "%m-%d-%y %H:%M",
+            },
+        },
+        "handlers": {
+            "console": {
+                "class": "logging.StreamHandler",
+                "level": "INFO",
+                "formatter": "standard",
+            },
+            "file": {
+                "class": "logging.FileHandler",
+                "filename": logfile,
+                "level": "INFO",
+                "mode": "w",
+                "formatter": "detailed",
+            }
+        },
+        "root": {
+            "level": "INFO",
+            "handlers": ["console", "file"]
+        },
+    })
+    logger = logging.getLogger(__name__)
+
+    kwargs = dict()
+    universe = mda.Universe(topology, trajectory)
+
+    if mass:
+        logger.info("Setting all bead masses to 1.0.")
+        universe.atoms.mass = 1.0
+
+    with mda.Writer(native_str(outfile), **kwargs) as rtf:
+        logger.info("Writing {}...".format(outfile))
+        rtf.write(universe, decl=not decl)

src/fluctmatch/fluctmatch/charmmfluctmatch.py

@@ -218,8 +218,9 @@ def initialize(self, restart=False):
 
             # Create and write initial internal coordinate files.
             logger.info("Determining the average bond distances...")
-            avg_bonds, std_bonds = fmutils.BondStats(
-                universe, func="both").run().result
+            avg_bonds = fmutils.BondAverage(universe).run().result
+            logger.info("Determining the bond distance fluctuations...")
+            std_bonds = fmutils.BondStd(universe, avg_bonds).run().result
             with mda.Writer(self.filenames["init_avg_ic"],
                             **self.kwargs) as table:
                 logger.info("Writing {}...".format(