Merge pull request #1 from Acellera/master

updating my fork

doc/source/index.rst

@@ -27,7 +27,7 @@ The user guide is a good place to start playing around.
 
 .. toctree::
     :maxdepth: 3
-
+    Installation <https://software.acellera.com/install-htmd.html>
     userguide
     api
     developers

htmd/__init__.py

@@ -4,6 +4,7 @@
 # No redistribution in whole or part
 #
 from htmd.version import version as _version
+from htmd.versionwarnings import _issueWarnings, _disableWarnings
 import __main__ as main
 import os.path
 from htmd.config import _config
@@ -37,3 +38,5 @@
             raise RuntimeError('Failed to execute the HTMD Config file {}.'.format(_config['configfile']))
         else:
             print('\nHTMD Config file {} executed.'.format(_config['configfile']))
+
+_issueWarnings()

htmd/adaptive/adaptive.py

@@ -198,7 +198,7 @@ def _writeInputs(self, simsframes, epoch=None):
 
         from htmd.parallelprogress import ParallelExecutor
         from htmd.config import _config
-        aprun = ParallelExecutor(n_jobs=_config['ncpus'])
+        aprun = ParallelExecutor(n_jobs=_config['njobs'])
         aprun(total=len(simsframes), desc='Writing inputs')(
             delayed(_writeInputsFunction)(i, f, epoch, self.inputpath, self.coorname) for i, f in enumerate(simsframes))
 

htmd/adaptive/adaptivegoal.py

@@ -68,7 +68,7 @@ class AdaptiveGoal(AdaptiveMD):
     skip : int, default=1
         Allows skipping of simulation frames to reduce data. i.e. skip=3 will only keep every third frame
     lag : int, default=1
-        The lagtime used to create the Markov model
+        The lagtime used to create the Markov model. Units are in frames.
     clustmethod : :class:`ClusterMixin <sklearn.base.ClusterMixin>` class, default=<class 'htmd.clustering.kcenters.KCenter'>
         Clustering algorithm used to cluster the contacts or distances
     method : str, default='1/Mc'

htmd/adaptive/adaptivegoaleg.py

@@ -68,7 +68,7 @@ class AdaptiveGoalEG(AdaptiveGoal):
     skip : int, default=1
         Allows skipping of simulation frames to reduce data. i.e. skip=3 will only keep every third frame
     lag : int, default=1
-        The lagtime used to create the Markov model
+        The lagtime used to create the Markov model. Units are in frames.
     clustmethod : :class:`ClusterMixin <sklearn.base.ClusterMixin>` class, default=<class 'htmd.clustering.kcenters.KCenter'>
         Clustering algorithm used to cluster the contacts or distances
     method : str, default='1/Mc'

htmd/adaptive/adaptiverun.py

@@ -70,7 +70,7 @@ class AdaptiveMD(AdaptiveBase):
     skip : int, default=1
         Allows skipping of simulation frames to reduce data. i.e. skip=3 will only keep every third frame
     lag : int, default=1
-        The lagtime used to create the Markov model
+        The lagtime used to create the Markov model. The units are in frames.
     clustmethod : :class:`ClusterMixin <sklearn.base.ClusterMixin>` class, default=<class 'htmd.clustering.kcenters.KCenter'>
         Clustering algorithm used to cluster the contacts or distances
     method : str, default='1/Mc'

htmd/adaptive/util.py

@@ -0,0 +1,66 @@
+import numpy as np
+import os
+import re
+
+parentregex = re.compile('e\d+s\d+_(e\d+s\d+)p(\d+)f(\d+)')
+epochregex = re.compile('^e(\d+)s')
+
+
+def getEpochTrajectoryDictionary(simlist):
+    from collections import defaultdict
+
+    simepochs = defaultdict(list)
+    for i, sim in enumerate(simlist):
+        name = os.path.basename(os.path.dirname(os.path.abspath(sim.trajectory[0])))
+        e = int(epochregex.findall(name)[0])
+        simepochs[e].append(i)
+
+    return simepochs
+
+def getEpochSimIdx(data, epoch):
+    idx = []
+    for i, sim in enumerate(data.simlist):
+        name = os.path.basename(os.path.dirname(os.path.abspath(sim.trajectory[0])))
+        e = int(epochregex.findall(name)[0])
+        if e == epoch:
+            idx.append(i)
+
+    return np.array(idx)
+
+def getParentSimIdxFrame(data, trajidx):
+    name = os.path.basename(os.path.dirname(os.path.abspath(data.simlist[trajidx].trajectory[0])))
+    res = parentregex.findall(name)
+    if len(res) == 0:
+        print('Failed to find parent simulation of {}. Assuming it\'s the first epoch. Will use the first frame as the state'.format(name))
+        return trajidx, 0
+
+    res = res[0]
+    parentname, parentpiece, parentframe = res
+    parentpiece = int(parentpiece)
+    parentframe = int(parentframe)
+
+    parentidx = None
+    for i, sim in enumerate(data.simlist):
+        name = os.path.basename(os.path.dirname(sim.trajectory[0]))
+        if name.startswith(parentname+'_'):
+            if parentidx is None:
+                parentidx = i
+            else:
+                raise RuntimeError('More than one simulation matches parent name!')
+
+    pieceframe = np.array([parentpiece, parentframe])
+    frameidx = np.where(np.all(data.trajectories[parentidx].reference == pieceframe, axis=1))[0]
+    if len(frameidx) != 1:
+        raise RuntimeError('Only one frame should match simpiece/frame combo')
+    frameidx = frameidx[0]
+
+    return parentidx, frameidx
+
+def updatingMean(oldmean, oldcount, newdata):
+    if oldcount == 0:
+        assert oldmean == 0
+        return np.mean(newdata)
+    newcount = oldcount + len(newdata)
+    return (oldmean + np.sum(newdata) / oldcount) * (oldcount / newcount)
+
+

htmd/builder/amber.py

@@ -65,7 +65,7 @@ def defaultAmberHome(teleap=None):
 def htmdAmberHome():
     """ Returns the location of the AMBER files distributed with HTMD"""
 
-    return os.path.abspath(os.path.join(home(), 'builder', 'amberfiles'))
+    return os.path.abspath(os.path.join(home(shareDir=True), 'builder', 'amberfiles'))
 
 
 def listFiles():
@@ -646,7 +646,7 @@ def _charmmLipid2Amber(mol):
         A new Molecule object with the membrane converted to AMBER
     """
 
-    resdict = _readcsvdict(os.path.join(home(), 'builder', 'charmmlipid2amber.csv'))
+    resdict = _readcsvdict(os.path.join(home(shareDir=True), 'builder', 'charmmlipid2amber.csv'))
 
     natoms = mol.numAtoms
     neworder = np.array(list(range(natoms)))  # After renaming the atoms and residues I have to reorder them
@@ -799,7 +799,7 @@ def _logParser(fname):
     return errors
 
 
-class TestAmberBuild(unittest.TestCase):
+class _TestAmberBuild(unittest.TestCase):
     currentResult = None  # holds last result object passed to run method
 
     def setUp(self):
@@ -867,7 +867,7 @@ def testWithProteinPrepare(self):
             _ = build(smol, ff=ffs, outdir=tmpdir)
 
             refdir = home(dataDir=join('test-amber-build', 'pp', pid))
-            TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)
+            _TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)
 
     def testWithoutProteinPrepare(self):
         from htmd.builder.solvate import solvate
@@ -884,7 +884,7 @@ def testWithoutProteinPrepare(self):
             _ = build(smol, ff=ffs, outdir=tmpdir)
 
             refdir = home(dataDir=join('test-amber-build', 'nopp', pid))
-            TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)
+            _TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)
 
     def testProteinLigand(self):
         from htmd.builder.solvate import solvate
@@ -913,7 +913,7 @@ def testProteinLigand(self):
         _ = build(smol, outdir=tmpdir, param=params, ionize=False)
 
         refdir = home(dataDir=join('test-amber-build', 'protLig', 'results'))
-        TestAmberBuild._compareResultFolders(refdir, tmpdir, '3PTB')
+        _TestAmberBuild._compareResultFolders(refdir, tmpdir, '3PTB')
 
         # # Test protein-ligand building
         # folder = home(dataDir='building-protein-ligand')
@@ -960,14 +960,14 @@ def test_customDisulfideBonds(self):
             _ = build(smol, ff=ffs, outdir=tmpdir, disulfide=disu)
 
             refdir = home(dataDir=join('test-amber-build', 'nopp', pid))
-            TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)
+            _TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)
 
             np.random.seed(1)
             tmpdir = os.path.join(self.testDir, 'withoutProtPrep', pid)
             _ = build(smol, ff=ffs, outdir=tmpdir)
 
             refdir = home(dataDir=join('test-amber-build', 'nopp', pid))
-            TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)
+            _TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)
 
     def test_customTeLeapImports(self):
         from htmd.builder.solvate import solvate
@@ -989,7 +989,7 @@ def test_customTeLeapImports(self):
             _ = build(smol, ff=ffs, outdir=tmpdir, teleapimports=teleapimports)
 
             refdir = home(dataDir=join('test-amber-build', 'nopp', pid))
-            TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)
+            _TestAmberBuild._compareResultFolders(refdir, tmpdir, pid)
 
 
 if __name__ == '__main__':

htmd/builder/charmm.py

@@ -28,6 +28,10 @@
 logger = logging.getLogger(__name__)
 
 
+def htmdCharmmHome():
+    """ Returns the location of the CHARMM files distributed with HTMD"""
+    return os.path.abspath(os.path.join(home(shareDir=True), 'builder', 'charmmfiles', ''))
+
 def listFiles():
     """ Lists all available Charmm topologies and parameter files
 
@@ -39,7 +43,7 @@ def listFiles():
 
     """
     from natsort import natsorted
-    charmmdir = path.join(home(), 'builder', 'charmmfiles', '')  # maybe just lookup current module?
+    charmmdir = htmdCharmmHome()
     topos = natsorted(glob(path.join(charmmdir, 'top', '*.rtf')))
     params = natsorted(glob(path.join(charmmdir, 'par', '*.prm')))
     streams = natsorted(glob(path.join(charmmdir, 'str', '*', '*.str')))
@@ -71,8 +75,7 @@ def search(key, name):
     --------
     >>> charmm.search(key='RESI', name = 'CHL1')  # doctest: +SKIP
     """
-    charmmdir = path.join(home(), 'builder', 'charmmfiles', '')
-    os.system('find {} -type f -exec grep -n "{} {}" {{}} +'.format(charmmdir, key, name))
+    os.system('find {} -type f -exec grep -n "{} {}" {{}} +'.format(htmdCharmmHome(), key, name))
 
 
 def defaultTopo():
@@ -82,7 +85,7 @@ def defaultTopo():
 
 def defaultParam():
     """ Returns the default parameters used by charmm.build """
-    return ['par/par_all36_prot_mod.prm', 'par/par_all36_lipid.prm', 'par/par_water_ions.prm', 'par/par_all36_cgenff.prm']
+    return ['par/par_all36_prot.prm', 'par/par_all36_lipid.prm', 'par/par_water_ions.prm', 'par/par_all36_cgenff.prm']
 
 
 def defaultStream():
@@ -91,7 +94,8 @@ def defaultStream():
 
 
 def build(mol, topo=None, param=None, stream=None, prefix='structure', outdir='./build', caps=None, ionize=True, saltconc=0,
-          saltanion=None, saltcation=None, disulfide=None, patches=None, noregen=None, aliasresidues=None, psfgen=None, execute=True, _clean=True):
+          saltanion=None, saltcation=None, disulfide=None, regenerate=['angles', 'dihedrals'], patches=None, noregen=None, 
+          aliasresidues=None, psfgen=None, execute=True, _clean=True):
     """ Builds a system for CHARMM
 
     Uses VMD and psfgen to build a system for CHARMM. Additionally it allows for ionization and adding of disulfide bridges.
@@ -107,7 +111,7 @@ def build(mol, topo=None, param=None, stream=None, prefix='structure', outdir='.
     param : list of str
         A list of parameter `prm` files.
         Use :func:`charmm.listFiles <htmd.builder.charmm.listFiles>` to get a list of available parameter files.
-        Default: ['par/par_all36_prot_mod.prm', 'par/par_all36_lipid.prm', 'par/par_water_ions.prm']
+        Default: ['par/par_all36_prot.prm', 'par/par_all36_lipid.prm', 'par/par_water_ions.prm']
     stream : list of str
         A list of stream `str` files containing topologies and parameters.
         Use :func:`charmm.listFiles <htmd.builder.charmm.listFiles>` to get a list of available stream files.
@@ -132,6 +136,9 @@ def build(mol, topo=None, param=None, stream=None, prefix='structure', outdir='.
     disulfide : list of pairs of atomselection strings
         If None it will guess disulfide bonds. Otherwise provide a list pairs of atomselection strings for each pair of
         residues forming the disulfide bridge.
+    regenerate : None or list of strings of: ['angles', 'dihedrals']
+        Disable angle/dihedral regeneration with `regenerate=None`, or enable it with `regenerate=['angles', 'diheldrals']` 
+        or just one of the two options with `regenerate=['angles']` or `regenerate=['diheldrals']`.
     patches : list of str
         Any further patches the user wants to apply
     noregen : list of str
@@ -159,7 +166,7 @@ def build(mol, topo=None, param=None, stream=None, prefix='structure', outdir='.
                  B: [serial 298 resid 58 resname CYS chain A segid 0]...
     >>> # More complex example
     >>> topos  = ['top/top_all36_prot.rtf', './benzamidine.rtf', 'top/top_water_ions.rtf']
-    >>> params = ['par/par_all36_prot_mod.prm', './benzamidine.prm', 'par/par_water_ions.prm']
+    >>> params = ['par/par_all36_prot.prm', './benzamidine.prm', 'par/par_water_ions.prm']
     >>> disu = [['segid P and resid 157', 'segid P and resid 13'], ['segid K and resid 1', 'segid K and resid 25']]
     >>> ar = {'SAPI24': 'SP24'}  # Alias large resnames to a short-hand version
     >>> molbuilt = charmm.build(mol, topo=topos, param=params, outdir='/tmp/build', saltconc=0.15, disulfide=disu, aliasresidues=ar)  # doctest: +SKIP
@@ -194,7 +201,7 @@ def build(mol, topo=None, param=None, stream=None, prefix='structure', outdir='.
     allparam = param.copy()
 
     # Splitting the stream files and adding them to the list of parameter and topology files
-    charmmdir = path.join(home(), 'builder', 'charmmfiles')
+    charmmdir = htmdCharmmHome()
     for s in stream:
         if s[0] != '.' and path.isfile(path.join(charmmdir, s)):
             s = path.join(charmmdir, s)
@@ -293,8 +300,9 @@ def build(mol, topo=None, param=None, stream=None, prefix='structure', outdir='.
         f.write('\n')
 
     # Regenerate angles and dihedrals
-    f.write('regenerate angles dihedrals\n')
-    f.write('\n')
+    if regenerate is not None:
+        f.write('regenerate {}\n'.format(' '.join(regenerate)))
+        f.write('\n')
 
     # Printing non-regenerable patches
     if len(noregenpatches) != 0:
@@ -347,7 +355,7 @@ def build(mol, topo=None, param=None, stream=None, prefix='structure', outdir='.
             newmol = ionizePlace(mol, anion, cation, anionatom, cationatom, nanion, ncation)
             # Redo the whole build but now with ions included
             return build(newmol, topo=alltopo, param=allparam, stream=[], prefix=prefix, outdir=outdir, ionize=False,
-                         caps=caps, execute=execute, saltconc=saltconc, disulfide=disulfide, patches=patches,
+                         caps=caps, execute=execute, saltconc=saltconc, disulfide=disulfide, regenerate=regenerate, patches=patches,
                          noregen=noregen, aliasresidues=aliasresidues, psfgen=psfgen, _clean=False)
     _checkFailedAtoms(molbuilt)
     _recoverProtonations(molbuilt)
@@ -652,7 +660,7 @@ def combine(prmlist, outfile):
     # Process parameter files
     prm_list = ["!COMMENTS\n", "!ATOMS\n", "BONDS\n", "ANGLES\n", "DIHEDRALS\n", "IMPROPER\n", "CMAP\n", "NONBONDED\n", "NBFIX\n", "HBOND\n"]
 
-    charmmdir = path.join(home(), 'builder', 'charmmfiles')
+    charmmdir = htmdCharmmHome()
     for myfile in prmlist:
         if myfile[0] != '.' and path.isfile(path.join(charmmdir, myfile)):
             myfile = path.join(charmmdir, myfile)
@@ -892,7 +900,7 @@ def _checkFailedAtoms(mol):
                         'Check log file for more details.'.format(idx))
 
 
-class TestCharmmBuild(TestCase):
+class _TestCharmmBuild(TestCase):
     def test_build(self):
         from moleculekit.molecule import Molecule
         from htmd.builder.solvate import solvate
@@ -916,7 +924,7 @@ def test_build(self):
                 np.random.seed(1)  # Needed for ions
                 smol = solvate(mol)
                 topos = ['top/top_all36_prot.rtf', 'top/top_water_ions.rtf']
-                params = ['par/par_all36_prot_mod.prm', 'par/par_water_ions.prm']
+                params = ['par/par_all36_prot.prm', 'par/par_water_ions.prm']
                 tmpdir = tempname()
                 _ = build(smol, topo=topos, param=params, outdir=tmpdir)
 
@@ -945,7 +953,7 @@ def test_customDisulfideBonds(self):
         np.random.seed(1)  # Needed for ions
         smol = solvate(mol)
         topos = ['top/top_all36_prot.rtf', 'top/top_water_ions.rtf']
-        params = ['par/par_all36_prot_mod.prm', 'par/par_water_ions.prm']
+        params = ['par/par_all36_prot.prm', 'par/par_water_ions.prm']
         disu = [['segid 1 and resid 110', 'segid 1 and resid 187'], ['segid 0 and resid 110', 'segid 0 and resid 187']]
         tmpdir = tempname()
         _ = build(smol, topo=topos, param=params, outdir=tmpdir, disulfide=disu)
@@ -985,7 +993,7 @@ def test_disulfideWithInsertion(self):
         np.random.seed(1)  # Needed for ions
         smol = solvate(mol)
         topos = ['top/top_all36_prot.rtf', 'top/top_water_ions.rtf']
-        params = ['par/par_all36_prot_mod.prm', 'par/par_water_ions.prm']
+        params = ['par/par_all36_prot.prm', 'par/par_water_ions.prm']
 
         smol.insertion[smol.resid == 42] = 'A'  # Adding an insertion to test that disulfide bonds with insertions work
         tmpdir = tempname()