Merge pull request #2 from orbeckst/master

packaged for standalone installation (Python 3 and Python 2)

INSTALL

@@ -0,0 +1,19 @@
+## Installation
+
+Clone repository or unpack the tar ball and install with
+[setuptools](http://pythonhosted.org/setuptools/index.html) (note: if
+you don't have setuptools installed you will need an internet
+connection so that the installation procedure can sownload the
+required files):
+
+    cd propka-3.1
+    python setup.py install --user
+
+This will install the `propka31` script in your executable directory,
+as configured for setuptools, for instance `~/.local/bin`. You can
+change the bin directory with the `--install-scripts` option. For
+example, in order to install in my `bin` directory in my home
+directory:
+
+    python setup.py install --user --install-scripts ~/bin
+

MANIFEST.in

@@ -0,0 +1,4 @@
+include README.md INSTALL 
+include ez_setup.py setup.py 
+include propka.cfg
+

README.md

@@ -1,45 +1,69 @@
 # PROPKA 3.1
 
-PROPKA predicts the pKa values of ionizable groups in
+PROPKA predicts the pKa values of ionizable groups in proteins
-proteins (version 3.0) and
+(version 3.0) and protein-ligand complexes (version 3.1)
-protein-ligand complexes (version 3.1)
+based on the 3D structure.
-based in the 3D structure.
 
 For proteins without ligands both version should produce the same result.
 
 The method is described in the following papers, which you should cite
 in publications:
 
-* Søndergaard, Chresten R., Mats HM Olsson, Michal Rostkowski, and Jan H. Jensen. "Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values." Journal of Chemical Theory and Computation 7, no. 7 (2011): 2284-2295.
+* Sondergaard, Chresten R., Mats HM Olsson, Michal Rostkowski, and Jan H. Jensen. "Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values." Journal of Chemical Theory and Computation 7, no. 7 (2011): 2284-2295.
 
-* Olsson, Mats HM, Chresten R. Søndergaard, Michal Rostkowski, and Jan H. Jensen. "PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions." Journal of Chemical Theory and Computation 7, no. 2 (2011): 525-537.
+* Olsson, Mats HM, Chresten R. Sondergaard, Michal Rostkowski, and Jan H. Jensen. "PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions." Journal of Chemical Theory and Computation 7, no. 2 (2011): 525-537.
 
 See [propka.ki.ku.dk](http://propka.ki.ku.dk/) for the PROPKA web server,
 using the [tutorial](http://propka.ki.ku.dk/~luca/wiki/index.php/PROPKA_3.1_Tutorial).
 
+## Modifications 
+
+This release of PROPKA 3.1 was modified by Oliver Beckstein
+<oliver.beckstein@asu.edu> from the released version.
+
+* Included patches from
+  https://github.com/schrodinger/propka-3.1/tree/python27-compat to
+  make it compatible with Python 2.7
+
+* Packaged for installation with setuptools.
+
+
 ## Installation
 
-No installation needed. Just clone and run.
+Clone repository or unpack the tar ball and install with
+[setuptools](http://pythonhosted.org/setuptools/index.html) (note: if
+you don't have setuptools installed you will need an internet
+connection so that the installation procedure can download the
+required files):
+
+    cd propka-3.1
+    python setup.py install --user
+
+This will install the `propka31` script in your executable directory,
+as configured for setuptools, for instance `~/.local/bin`. You can
+change the bin directory with the `--install-scripts` option. For
+example, in order to install in my `bin` directory in my home
+directory:
+
+    python setup.py install --user --install-scripts ~/bin
+
 
 ## Requirements
 
-* Python 3.1 or higher
+* Python 2.7 or higher or Python 3.1 or higher 
 
 ## Getting started
 
 1. Clone the code from GitHub
-2. Run 'propka.py' with a .pdb file (see Examples)
+2. `python setup.py install --user`
+2. Run `propka31` with a .pdb file (see Examples)
 
 ## Examples
 
 Calculate using pdb file
 
-    ./propka.py 1hpx.pdb
+    propka31 1hpx.pdb
 
-If for some reason your setup with python3.1+ is
-not located in '/usr/bin/python3', run the script
-
-    python3.2 propka.py 1hpx.pdb
 
 ## Testing (for developers)
 
@@ -49,9 +73,9 @@ Please run `Tests/runtest.py/` after changes before pushing commits.
 
 Please cite these references in publications:
 
-* Søndergaard, Chresten R., Mats HM Olsson, Michal Rostkowski, and Jan H. Jensen. "Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values." Journal of Chemical Theory and Computation 7, no. 7 (2011): 2284-2295.
+* Sondergaard, Chresten R., Mats HM Olsson, Michal Rostkowski, and Jan H. Jensen. "Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values." Journal of Chemical Theory and Computation 7, no. 7 (2011): 2284-2295.
 
-* Olsson, Mats HM, Chresten R. Søndergaard, Michal Rostkowski, and Jan H. Jensen. "PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions." Journal of Chemical Theory and Computation 7, no. 2 (2011): 525-537.
+* Olsson, Mats HM, Chresten R. Sondergaard, Michal Rostkowski, and Jan H. Jensen. "PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions." Journal of Chemical Theory and Computation 7, no. 2 (2011): 525-537.
 
 
 

Tests/runtest.py

@@ -1,9 +1,13 @@
-#!/usr/bin/python3
+#!/usr/bin/env python
 
 """ Run test for test pdbs """
 
+from __future__ import division
+from __future__ import print_function
+
 from subprocess import call
 import os, re
+import sys
 
 pdbs = ['1FTJ-Chain-A',
         '1HPX',
@@ -14,7 +18,7 @@ for pdb in pdbs:
   print('RUNNING '+pdb)
 
   # Run pka calculation
-  call(['../propka.py','pdb/'+pdb+'.pdb'], stdout = open(pdb+'.out', 'w+'))
+  call([sys.executable, '../scripts/propka31.py','pdb/'+pdb+'.pdb'], stdout = open(pdb+'.out', 'w+'))
 
   # Test pka predictiona
   result = open('results/'+pdb+'.dat','r')

ez_setup.py

@@ -0,0 +1,258 @@
+#!python
+"""Bootstrap setuptools installation
+
+If you want to use setuptools in your package's setup.py, just include this
+file in the same directory with it, and add this to the top of your setup.py::
+
+    from ez_setup import use_setuptools
+    use_setuptools()
+
+If you want to require a specific version of setuptools, set a download
+mirror, or use an alternate download directory, you can do so by supplying
+the appropriate options to ``use_setuptools()``.
+
+This file can also be run as a script to install or upgrade setuptools.
+"""
+import os
+import shutil
+import sys
+import tempfile
+import tarfile
+import optparse
+import subprocess
+
+from distutils import log
+
+try:
+    from site import USER_SITE
+except ImportError:
+    USER_SITE = None
+
+DEFAULT_VERSION = "0.9.6"
+DEFAULT_URL = "https://pypi.python.org/packages/source/s/setuptools/"
+
+def _python_cmd(*args):
+    args = (sys.executable,) + args
+    return subprocess.call(args) == 0
+
+def _install(tarball, install_args=()):
+    # extracting the tarball
+    tmpdir = tempfile.mkdtemp()
+    log.warn('Extracting in %s', tmpdir)
+    old_wd = os.getcwd()
+    try:
+        os.chdir(tmpdir)
+        tar = tarfile.open(tarball)
+        _extractall(tar)
+        tar.close()
+
+        # going in the directory
+        subdir = os.path.join(tmpdir, os.listdir(tmpdir)[0])
+        os.chdir(subdir)
+        log.warn('Now working in %s', subdir)
+
+        # installing
+        log.warn('Installing Setuptools')
+        if not _python_cmd('setup.py', 'install', *install_args):
+            log.warn('Something went wrong during the installation.')
+            log.warn('See the error message above.')
+            # exitcode will be 2
+            return 2
+    finally:
+        os.chdir(old_wd)
+        shutil.rmtree(tmpdir)
+
+
+def _build_egg(egg, tarball, to_dir):
+    # extracting the tarball
+    tmpdir = tempfile.mkdtemp()
+    log.warn('Extracting in %s', tmpdir)
+    old_wd = os.getcwd()
+    try:
+        os.chdir(tmpdir)
+        tar = tarfile.open(tarball)
+        _extractall(tar)
+        tar.close()
+
+        # going in the directory
+        subdir = os.path.join(tmpdir, os.listdir(tmpdir)[0])
+        os.chdir(subdir)
+        log.warn('Now working in %s', subdir)
+
+        # building an egg
+        log.warn('Building a Setuptools egg in %s', to_dir)
+        _python_cmd('setup.py', '-q', 'bdist_egg', '--dist-dir', to_dir)
+
+    finally:
+        os.chdir(old_wd)
+        shutil.rmtree(tmpdir)
+    # returning the result
+    log.warn(egg)
+    if not os.path.exists(egg):
+        raise IOError('Could not build the egg.')
+
+
+def _do_download(version, download_base, to_dir, download_delay):
+    egg = os.path.join(to_dir, 'setuptools-%s-py%d.%d.egg'
+                       % (version, sys.version_info[0], sys.version_info[1]))
+    if not os.path.exists(egg):
+        tarball = download_setuptools(version, download_base,
+                                      to_dir, download_delay)
+        _build_egg(egg, tarball, to_dir)
+    sys.path.insert(0, egg)
+    import setuptools
+    setuptools.bootstrap_install_from = egg
+
+
+def use_setuptools(version=DEFAULT_VERSION, download_base=DEFAULT_URL,
+                   to_dir=os.curdir, download_delay=15):
+    # making sure we use the absolute path
+    to_dir = os.path.abspath(to_dir)
+    was_imported = 'pkg_resources' in sys.modules or \
+        'setuptools' in sys.modules
+    try:
+        import pkg_resources
+    except ImportError:
+        return _do_download(version, download_base, to_dir, download_delay)
+    try:
+        pkg_resources.require("setuptools>=" + version)
+        return
+    except pkg_resources.VersionConflict:
+        e = sys.exc_info()[1]
+        if was_imported:
+            sys.stderr.write(
+            "The required version of setuptools (>=%s) is not available,\n"
+            "and can't be installed while this script is running. Please\n"
+            "install a more recent version first, using\n"
+            "'easy_install -U setuptools'."
+            "\n\n(Currently using %r)\n" % (version, e.args[0]))
+            sys.exit(2)
+        else:
+            del pkg_resources, sys.modules['pkg_resources']    # reload ok
+            return _do_download(version, download_base, to_dir,
+                                download_delay)
+    except pkg_resources.DistributionNotFound:
+        return _do_download(version, download_base, to_dir,
+                            download_delay)
+
+
+def download_setuptools(version=DEFAULT_VERSION, download_base=DEFAULT_URL,
+                        to_dir=os.curdir, delay=15):
+    """Download setuptools from a specified location and return its filename
+
+    `version` should be a valid setuptools version number that is available
+    as an egg for download under the `download_base` URL (which should end
+    with a '/'). `to_dir` is the directory where the egg will be downloaded.
+    `delay` is the number of seconds to pause before an actual download
+    attempt.
+    """
+    # making sure we use the absolute path
+    to_dir = os.path.abspath(to_dir)
+    try:
+        from urllib.request import urlopen
+    except ImportError:
+        from urllib2 import urlopen
+    tgz_name = "setuptools-%s.tar.gz" % version
+    url = download_base + tgz_name
+    saveto = os.path.join(to_dir, tgz_name)
+    src = dst = None
+    if not os.path.exists(saveto):  # Avoid repeated downloads
+        try:
+            log.warn("Downloading %s", url)
+            src = urlopen(url)
+            # Read/write all in one block, so we don't create a corrupt file
+            # if the download is interrupted.
+            data = src.read()
+            dst = open(saveto, "wb")
+            dst.write(data)
+        finally:
+            if src:
+                src.close()
+            if dst:
+                dst.close()
+    return os.path.realpath(saveto)
+
+
+def _extractall(self, path=".", members=None):
+    """Extract all members from the archive to the current working
+       directory and set owner, modification time and permissions on
+       directories afterwards. `path' specifies a different directory
+       to extract to. `members' is optional and must be a subset of the
+       list returned by getmembers().
+    """
+    import copy
+    import operator
+    from tarfile import ExtractError
+    directories = []
+
+    if members is None:
+        members = self
+
+    for tarinfo in members:
+        if tarinfo.isdir():
+            # Extract directories with a safe mode.
+            directories.append(tarinfo)
+            tarinfo = copy.copy(tarinfo)
+            tarinfo.mode = 448  # decimal for oct 0700
+        self.extract(tarinfo, path)
+
+    # Reverse sort directories.
+    if sys.version_info < (2, 4):
+        def sorter(dir1, dir2):
+            return cmp(dir1.name, dir2.name)
+        directories.sort(sorter)
+        directories.reverse()
+    else:
+        directories.sort(key=operator.attrgetter('name'), reverse=True)
+
+    # Set correct owner, mtime and filemode on directories.
+    for tarinfo in directories:
+        dirpath = os.path.join(path, tarinfo.name)
+        try:
+            self.chown(tarinfo, dirpath)
+            self.utime(tarinfo, dirpath)
+            self.chmod(tarinfo, dirpath)
+        except ExtractError:
+            e = sys.exc_info()[1]
+            if self.errorlevel > 1:
+                raise
+            else:
+                self._dbg(1, "tarfile: %s" % e)
+
+
+def _build_install_args(options):
+    """
+    Build the arguments to 'python setup.py install' on the setuptools package
+    """
+    install_args = []
+    if options.user_install:
+        if sys.version_info < (2, 6):
+            log.warn("--user requires Python 2.6 or later")
+            raise SystemExit(1)
+        install_args.append('--user')
+    return install_args
+
+def _parse_args():
+    """
+    Parse the command line for options
+    """
+    parser = optparse.OptionParser()
+    parser.add_option(
+        '--user', dest='user_install', action='store_true', default=False,
+        help='install in user site package (requires Python 2.6 or later)')
+    parser.add_option(
+        '--download-base', dest='download_base', metavar="URL",
+        default=DEFAULT_URL,
+        help='alternative URL from where to download the setuptools package')
+    options, args = parser.parse_args()
+    # positional arguments are ignored
+    return options
+
+def main(version=DEFAULT_VERSION):
+    """Install or upgrade setuptools and EasyInstall"""
+    options = _parse_args()
+    tarball = download_setuptools(download_base=options.download_base)
+    return _install(tarball, _build_install_args(options))
+
+if __name__ == '__main__':
+    sys.exit(main())

propka.py

@@ -1,23 +0,0 @@
-#!/usr/bin/python3
-
-import string,re,sys,os,math
-import Source.lib, Source.molecular_container 
-
-def main():
-    """
-    Reads in structure files, calculates pKa values, and prints pKa files
-    """
-    # loading options, flaggs and arguments
-    options, pdbfiles = Source.lib.loadOptions()
-
-    for pdbfile in pdbfiles:
-        my_molecule = Source.molecular_container.Molecular_container(pdbfile, options)
-        my_molecule.calculate_pka()
-        my_molecule.write_pka()
-
-
-if __name__ == '__main__': 
-    #import cProfile
-    #cProfile.run('main()',sort=1)
-    main()
-

propka/__init__.py

@@ -1,3 +1,5 @@
-#
+# propka 3.1
+# (Module renamed from original Source to propka to facilitate setuptools installation)
+
 __all__ = ['coupled_residues', 'lib', 'parameters', 'residue', 'bonds', 'debug', 'ligand', 'pdb', 'calculator', 'determinants', 'mutate', 'protein', 'chain', 'iterative', 'output', 'protonate']
 

propka/atom.py

@@ -1,4 +1,8 @@
-import string, Source.lib, Source.group
+
+from __future__ import division
+from __future__ import print_function
+
+import string, propka.lib, propka.group
 
 
 class Atom:
@@ -8,17 +12,17 @@ class Atom:
 
     def __init__(self, line=None, verbose=False):
 
-        self.set_properties(line)     
+        self.set_properties(line)
 
-        self.residue_label = "%-3s%4d%2s" % (self.name,self.resNumb, self.chainID)        
+        self.residue_label = "%-3s%4d%2s" % (self.name,self.resNumb, self.chainID)
 
         self.groups_extracted = 0
         self.group = None
         self.group_type = None
         self.number_of_bonded_elements = {}
-        
+
         self.cysteine_bridge = False
-        
+
         self.bonded_atoms = []
         self.residue = None
         self.conformation_container = None
@@ -39,7 +43,7 @@ class Atom:
         self.sybyl_type = ''
         self.sybyl_assigned = False
         self.marvin_pka = False
-        
+
 
         return
 
@@ -122,13 +126,13 @@ class Atom:
                     return True
 
         return False
-        
+
 
 
     def setProperty(self,
-                    numb    = None, 
+                    numb    = None,
-                    name    = None, 
+                    name    = None,
-                    resName = None, 
+                    resName = None,
                     chainID = None,
                     resNumb = None,
                     x       = None,
@@ -157,7 +161,7 @@ class Atom:
         """
         making a copy of this atom
         """
-        
+
         newAtom = Atom()
         newAtom.type = self.type
         newAtom.numb = self.numb
@@ -191,11 +195,11 @@ class Atom:
 
             if self.group.titratable:
                 model_pka = '%6.2f'%self.group.model_pka
-       
+
 
         str  = "%-6s%5d %s %s%2s%4d%12.3lf%8.3lf%8.3lf%6s%6s \n" % (self.type.upper(),
                                                                     self.numb,
-                                                                    Source.lib.makeTidyAtomLabel(self.name,
+                                                                    propka.lib.makeTidyAtomLabel(self.name,
                                                                                                  self.element),
                                                                     self.resName,
                                                                     self.chainID,
@@ -213,7 +217,7 @@ class Atom:
 
     def make_conect_line(self):
         res = 'CONECT%5d'%self.numb
-        
+
         # extract and sort numbers of bonded residues
         bonded = []
         for atom in self.bonded_atoms:
@@ -226,9 +230,9 @@ class Atom:
         res += '\n'
         return res
 
-    
+
     def get_input_parameters(self):
-        """ Method for getting the input parameters stored in the 
+        """ Method for getting the input parameters stored in the
         occupancy and b-factor fields in input files"""
 
         # Set the group type
@@ -258,7 +262,7 @@ class Atom:
 
             # try to initialise the group
             try:
-                exec('self.group = Source.group.%s_group(self)'%self.occ)
+                exec('self.group = propka.group.%s_group(self)'%self.occ)
             except:
                 raise Exception('%s in input_file is not recognized as a group'%self.occ)
 
@@ -270,7 +274,7 @@ class Atom:
         # set occ and beta to standard values
         self.occ = '1.00'
         self.beta = '0.00'
-        
+
 
 
         return
@@ -282,7 +286,7 @@ class Atom:
         # making string
         str  = "%-6s%5d %s %s%2s%4d%12.3lf%8.3lf%8.3lf%6s%6s\n" % (self.type.upper(),
                                                                    self.numb,
-                                                                   Source.lib.makeTidyAtomLabel(self.name,
+                                                                   propka.lib.makeTidyAtomLabel(self.name,
                                                                                                 self.element),
                                                                    self.resName,
                                                                    self.chainID,
@@ -292,17 +296,17 @@ class Atom:
                                                                    self.z,
                                                                    self.occ,
                                                                    self.beta)
-        
+
 
         return str
 
 
 
     def make_mol2_line(self,id):
-#1	S1     3.6147     2.0531     1.4795	S.3	1	noname	-0.1785
+#1      S1     3.6147     2.0531     1.4795     S.3     1       noname  -0.1785
         # making string
         str  = "%-4d %-4s %10.4f %10.4f %10.4f %6s %6d %10s %10.4f\n" % (id,
-                                                                         Source.lib.makeTidyAtomLabel(self.name,
+                                                                         propka.lib.makeTidyAtomLabel(self.name,
                                                                                                       self.element),
                                                                          self.x,
                                                                          self.y,
@@ -311,7 +315,7 @@ class Atom:
                                                                          self.resNumb,
                                                                          self.resName,
                                                                          0.0)#self.charge)
-        
+
 
         return str
 
@@ -346,7 +350,7 @@ class Atom:
         # making string
         str  = "ATOM "
         str += "%6d" % (numb)
-        str += " %s" % (Source.lib.makeTidyAtomLabel(name,self.element))
+        str += " %s" % (propka.lib.makeTidyAtomLabel(name,self.element))
         str += " %s" % (resName)
         str += "%2s" % (chainID)
         str += "%4d" % (resNumb)
@@ -363,25 +367,25 @@ class Atom:
         """
         Returns a 'tidier' atom label for printing the new pdbfile
         """
-        return Source.lib.makeTidyAtomLabel(self.name,self.element)
+        return propka.lib.makeTidyAtomLabel(self.name,self.element)
 
 
     def __str__(self):
         return '%5d-%4s %5d-%3s (%1s) [%8.3f %8.3f %8.3f] %s' %(self.numb, self.name, self.resNumb, self.resName, self.chainID, self.x, self.y, self.z,self.element)
-            
+
 
 #    def get_element(self):
 #        """ try to extract element if not already done"""
 #        if self.element == '':
 #            self.element = self.name.strip(string.digits)[0]
 #        return self.element
-    
+
 
     def set_residue(self, residue):
         """ Makes a references to the parent residue"""
         if self.residue == None:
             self.residue = residue
 
- 
+
 
 

propka/bonds.py

@@ -1,14 +1,19 @@
-import pickle,sys,os,math,Source.calculations
 
+from __future__ import division
+from __future__ import print_function
+
+import pickle,sys,os,math,propka.calculations
+
+import pkg_resources
 
 class bondmaker:
     def __init__(self):
-        
+
         # predefined bonding distances
         self.distances = {
             'S-S':2.5,
             'F-F':1.7}
-    
+
         self.distances_squared = {}
         for k in self.distances.keys():
             self.distances_squared[k]=self.distances[k]*self.distances[k]
@@ -18,12 +23,11 @@ class bondmaker:
 
         self.H_dist_squared  = self.H_dist * self.H_dist
         self.default_dist_squared = self.default_dist * self.default_dist
-        
+
-        self.max_sq_distance = max(list(self.distances_squared.values())+[self.default_dist_squared]) 
+        self.max_sq_distance = max(list(self.distances_squared.values())+[self.default_dist_squared])
-        
+
         # protein bonding data
-        path = os.path.split(__file__)[0]
+        self.data_file_name = pkg_resources.resource_filename(__name__, 'protein_bonds.dat')
-        self.data_file_name = os.path.join(path,'protein_bonds.dat')
 
         data = open(self.data_file_name,'rb')
         self.protein_bonds = pickle.load(data)
@@ -104,13 +108,13 @@ class bondmaker:
 
 
     def find_bonds_for_protein(self, protein):
-        """ Finds bonds proteins based on the way atoms 
+        """ Finds bonds proteins based on the way atoms
         normally bond in proteins"""
 
         print('++++ Side chains ++++')
         # side chains
         for chain in protein.chains:
-            for residue in chain.residues: 
+            for residue in chain.residues:
                 if residue.resName.replace(' ','') not in ['N+','C-']:
                     self.find_bonds_for_side_chain(residue.atoms)
 
@@ -145,7 +149,7 @@ class bondmaker:
             if atom1.name == 'SG':
                 for atom2 in cys2.atoms:
                     if atom2.name == 'SG':
-                        if Source.calculations.squared_distance(atom1,atom2) < self.SS_dist_squared:
+                        if propka.calculations.squared_distance(atom1,atom2) < self.SS_dist_squared:
                             self.make_bond(atom1, atom2)
 
 
@@ -174,7 +178,7 @@ class bondmaker:
             if atom1.name == 'C':
                 for atom2 in residue2.atoms:
                     if atom2.name == 'N':
-                        if Source.calculations.squared_distance(atom1,atom2) < self.default_dist_squared:
+                        if propka.calculations.squared_distance(atom1,atom2) < self.default_dist_squared:
                             self.make_bond(atom1, atom2)
 
         return
@@ -209,7 +213,7 @@ class bondmaker:
                     for atom2 in atoms:
                         if atom2.name in self.protein_bonds[atom1.resName][atom1.name]:
                             self.make_bond(atom1,atom2)
-                        
+
         return
 
 
@@ -270,7 +274,7 @@ class bondmaker:
 
                 if atoms[i] in atoms[j].bonded_atoms:
                     continue
-        
+
                 if self.check_distance(atoms[i], atoms[j]):
                     self.make_bond(atoms[i],atoms[j])
                     # di-sulphide bonds
@@ -283,14 +287,14 @@ class bondmaker:
 
 
     def check_distance(self, atom1, atom2):
-        sq_dist = Source.calculations.squared_distance(atom1, atom2)
+        sq_dist = propka.calculations.squared_distance(atom1, atom2)
-        
+
         if sq_dist > self.max_sq_distance:
             return False
-        
+
         key = '%s-%s'%(atom1.element,atom2.element)
         h_count = key.count('H')
-        
+
         if sq_dist <  self.H_dist_squared and h_count==1:
             return True
         if sq_dist < self.default_dist_squared and h_count==0:
@@ -299,8 +303,8 @@ class bondmaker:
             if sq_dist < self.distances_squared[key]:
                 return True
 
-        
+
- 
+
         return False
 
     def find_bonds_for_molecules_using_boxes(self, molecules):
@@ -320,7 +324,7 @@ class bondmaker:
 
     def find_bonds_for_atoms_using_boxes(self, atoms):
         """ Finds all bonds for a list of atoms"""
-        
+
         box_size = 2.5
 
         # find min and max coordinates
@@ -348,11 +352,13 @@ class bondmaker:
         #print('x range: [%6.2f;%6.2f] %6.2f'%(xmin,xmax,xlen))
         #print('y range: [%6.2f;%6.2f] %6.2f'%(ymin,ymax,ylen))
         #print('z range: [%6.2f;%6.2f] %6.2f'%(zmin,zmax,zlen))
-        
+
         # how many boxes do we need in each dimension?
-        self.no_box_x = max(1, math.ceil(xlen/box_size))
+        # NOTE: math.ceil() returns an int in python3 and a float in python2,
-        self.no_box_y = max(1, math.ceil(ylen/box_size))
+        # so we need to cast it to int for range() to work.
-        self.no_box_z = max(1, math.ceil(zlen/box_size))
+        self.no_box_x = max(1, int(math.ceil(xlen/box_size)))
+        self.no_box_y = max(1, int(math.ceil(ylen/box_size)))
+        self.no_box_z = max(1, int(math.ceil(zlen/box_size)))
 
         #print('No. box x: %6.2f'%self.no_box_x)
         #print('No. box y: %6.2f'%self.no_box_y)
@@ -372,12 +378,12 @@ class bondmaker:
             z = math.floor((atom.z-zmin)/box_size)
             self.put_atom_in_box(x,y,z,atom)
 
-        # assign bonds 
+        # assign bonds
         keys = self.boxes.keys()
         for key in keys:
             self.find_bonds_for_atoms(self.boxes[key])
 
-            
+
 
         return
 
@@ -386,7 +392,7 @@ class bondmaker:
         # one side of the x,y,z box in each dimension
         for bx in [x,x+1]:
             for by in [y,y+1]:
-                for bz in [z,z+1]:   
+                for bz in [z,z+1]:
                     key = self.box_key(bx,by,bz)
                     if key in self.boxes.keys():
                         self.boxes[key].append(atom)
@@ -397,13 +403,13 @@ class bondmaker:
 
     def box_key(self, x, y, z):
         return '%d-%d-%d'%(x,y,z)
-        
+
 
     def has_bond(self, atom1, atom2):
         if atom1 in atom2.bonded_atoms or atom2 in atom1.bonded_atoms:
             return True
         return False
-        
+
 
     def make_bond(self, atom1, atom2):
         """ Makes a bond between atom1 and atom2 """
@@ -413,13 +419,13 @@ class bondmaker:
         #print('making bond for',atom1,atom2)
         if not atom1 in atom2.bonded_atoms:
             atom2.bonded_atoms.append(atom1)
-        
+
         if not atom2 in atom1.bonded_atoms:
-            atom1.bonded_atoms.append(atom2)       
+            atom1.bonded_atoms.append(atom2)
 
         return
 
-        
+
     def generate_protein_bond_dictionary(self, atoms):
 
         for atom in atoms:
@@ -433,25 +439,25 @@ class bondmaker:
                         not name_j in self.backbone_atoms:
                     if not name_i in self.terminal_oxygen_names and\
                             not name_j in self.terminal_oxygen_names:
-                    
+
                         if not resi_i in list(self.protein_bonds.keys()):
                             self.protein_bonds[resi_i] = {}
                         if not name_i in self.protein_bonds[resi_i]:
                             self.protein_bonds[resi_i][name_i] = []
-                
+
                         if not name_j in self.protein_bonds[resi_i][name_i]:
                             self.protein_bonds[resi_i][name_i].append(name_j)
-            
+
 
                         if not resi_j in list(self.protein_bonds.keys()):
                             self.protein_bonds[resi_j] = {}
                         if not name_j in self.protein_bonds[resi_j]:
                             self.protein_bonds[resi_j][name_j] = []
-                            
+
                         if not name_i in self.protein_bonds[resi_j][name_j]:
                             self.protein_bonds[resi_j][name_j].append(name_i)
 
-        
+
         return
 
 
@@ -467,10 +473,10 @@ if __name__ == '__main__':
     if not os.path.isfile(filename):
         print('Error: Could not find \"%s\"'%filename)
         sys.exit(1)
-   
+
     pdblist = pdb.readPDB(filename)
     my_protein = protein.Protein(pdblist,'test.pdb')
-    
+
     for chain in my_protein.chains:
         for residue in chain.residues:
             residue.atoms = [atom for atom in residue.atoms if atom.element != 'H']

propka/calculations.py

@@ -1,6 +1,8 @@
 
+from __future__ import division
+from __future__ import print_function
 
-import math, Source.protonate, Source.bonds,copy, sys
+import math, propka.protonate, propka.bonds,copy, sys
 
 
 #
@@ -17,9 +19,9 @@ def setup_bonding_and_protonation(parameters, molecular_container):
     # apply information on pi electrons
     my_bond_maker.add_pi_electron_information(molecular_container)
 
-    # Protonate atoms 
+    # Protonate atoms
     if molecular_container.options.protonate_all:
-        my_protonator = Source.protonate.Protonate(verbose=False)
+        my_protonator = propka.protonate.Protonate(verbose=False)
         my_protonator.protonate(molecular_container)
 
 
@@ -29,7 +31,7 @@ def setup_bonding_and_protonation(parameters, molecular_container):
 
 def setup_bonding(parameters, molecular_container):
     # make bonds
-    my_bond_maker = Source.bonds.bondmaker()
+    my_bond_maker = propka.bonds.bondmaker()
     my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)
 
     return my_bond_maker
@@ -46,11 +48,11 @@ def set_ligand_atom_names(molecular_container):
 # The following methods imitates propka3.0 protonation!
 #
 def setup_bonding_and_protonation_30_style(parameters, molecular_container):
-    # Protonate atoms 
+    # Protonate atoms
     protonate_30_style(molecular_container)
 
     # make bonds
-    my_bond_maker = Source.bonds.bondmaker()
+    my_bond_maker = propka.bonds.bondmaker()
     my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)
 
     return
@@ -87,7 +89,7 @@ def protonate_30_style(molecular_container):
                     residue = []
             if atom.type=='atom':
                 residue.append(atom)
-                    
+
     return
 
 def addArgHydrogen(residue):
@@ -117,7 +119,7 @@ def addArgHydrogen(residue):
     H4.name = "HN3"
     H5 = protonateDirection([NH2, NE, H1])
     H5.name = "HN4"
-    
+
     return [H1,H2,H3,H4,H5]
 
 def addHisHydrogen(residue):
@@ -179,7 +181,7 @@ def addAmdHydrogen(residue):
             O = atom
         elif (atom.resName == "GLN" and atom.name == "NE2") or (atom.resName == "ASN" and atom.name == "ND2"):
             N = atom
-     
+
     if C == None or O == None or N == None:
         str = "Did not find N, C and/or O in %s%4d - in %s" % (atom.resName, atom.resNumb, "addAmdHydrogen()")
         print(str)
@@ -189,20 +191,20 @@ def addAmdHydrogen(residue):
     H1.name = "HN1"
     H2 = protonateAverageDirection([N, C, O])
     H2.name = "HN2"
-    
+
     return
 
 def addBackBoneHydrogen(residue, O, C):
     """
-    Adds hydrogen backbone atoms to residues according to the old way; dR is wrong for the N-terminus 
+    Adds hydrogen backbone atoms to residues according to the old way; dR is wrong for the N-terminus
     (i.e. first residue) but it doesn't affect anything at the moment. Could be improved, but works
     for now.
     """
-    
+
     new_C     = None
     new_O     = None
     N     = None
-    
+
 
     for atom in residue:
         if atom.name == "N":
@@ -211,20 +213,20 @@ def addBackBoneHydrogen(residue, O, C):
             new_C = atom
         if atom.name == "O":
             new_O = atom
-        
+
-        
+
-    
+
 
     if None in [C,O,N]:
         return [new_O,new_C]
 
-  
+
     if N.resName == "PRO":
         """ PRO doesn't have an H-atom; do nothing """
     else:
         H = protonateDirection([N, O, C])
         H.name = "H"
-        
+
     return [new_O,new_C]
 
 
@@ -247,7 +249,7 @@ def protonateDirection(list):
     y = X1.y + dY/length
     z = X1.z + dZ/length
 
-    
+
     H  = make_new_H(X1,x,y,z)
     H.name = "H"
 
@@ -277,7 +279,7 @@ def protonateAverageDirection(list):
     H  = make_new_H(X1,x,y,z)
     H.name = "H"
 
-    
+
 
     return H
 
@@ -290,7 +292,7 @@ def protonateSP2(list):
     X1 = list[0]
     X2 = list[1]
     X3 = list[2]
-    
+
     dX = (X1.x + X3.x)*0.5 - X2.x
     dY = (X1.y + X3.y)*0.5 - X2.y
     dZ = (X1.z + X3.z)*0.5 - X2.z
@@ -301,16 +303,16 @@ def protonateSP2(list):
 
     H  = make_new_H(X2,x,y,z)
     H.name = "H"
-    
+
     return H
 
 
 def make_new_H(atom, x,y,z):
 
-    new_H = Source.atom.Atom()
+    new_H = propka.atom.Atom()
-    new_H.setProperty(numb    = None, 
+    new_H.setProperty(numb    = None,
-                      name    = 'H%s'%atom.name[1:], 
+                      name    = 'H%s'%atom.name[1:],
-                      resName = atom.resName, 
+                      resName = atom.resName,
                       chainID = atom.chainID,
                       resNumb = atom.resNumb,
                       x       = x,
@@ -329,7 +331,7 @@ def make_new_H(atom, x,y,z):
     new_H.number_of_pi_electrons_in_double_and_triple_bonds = 0
 
     atom.bonded_atoms.append(new_H)
-    atom.conformation_container.add_atom(new_H)        
+    atom.conformation_container.add_atom(new_H)
 
     return new_H
 
@@ -356,7 +358,7 @@ def radial_volume_desolvation(parameters, group):
         # ignore atoms in the same residue
         if atom.resNumb == group.atom.resNumb and atom.chainID == group.atom.chainID:
             continue
-        
+
         sq_dist = squared_distance(group, atom)
 
         # desolvation
@@ -375,7 +377,7 @@ def radial_volume_desolvation(parameters, group):
         # buried
         if sq_dist < parameters.buried_cutoff_squared:
             group.Nmass += 1
-            
+
     group.buried = calculate_weight(parameters, group.Nmass)
     scale_factor = calculate_scale_factor(parameters, group.buried)
     volume_after_allowance = max(0.00, volume-parameters.desolvationAllowance)
@@ -385,9 +387,9 @@ def radial_volume_desolvation(parameters, group):
     # Elocl, Nlocl -> reorganisation energy (count backbone hydorgen bond acceptors, C=O)
 
 
-                
+
     #print('%s %5.2f %5.2f %4d'%(group, group.buried, group.Emass, group.Nmass))
-    return 
+    return
 
 
 
@@ -395,7 +397,7 @@ def contactDesolvation(parameters, group):
     """
     calculates the desolvation according to the Contact Model, the old default
     """
-    
+
     local_radius = {'ASP': 4.5,
                     'GLU': 4.5,
                     'HIS': 4.5,
@@ -405,7 +407,7 @@ def contactDesolvation(parameters, group):
                     'ARG': 5.0,
                     'C-': 4.5,
                     'N+': 4.5}
-    
+
     all_atoms = group.atom.conformation_container.get_non_hydrogen_atoms()
     if residue.resName in version.desolvationRadii:
         local_cutoff = version.desolvationRadii[residue.resName]
@@ -434,7 +436,7 @@ def contactDesolvation(parameters, group):
         # Note, there will be an unforseen problem: e.g. if one residue has Nmass > Nmax and
         # the other Nmass < Nmax, the Npair will not be Nmass1 + Nmass2!
         residue.buried = calculateWeight(residue.Nmass)
-    
+
         return 0.00, 0.00, 0.00, 0.00
 
 
@@ -455,10 +457,10 @@ def calculate_weight(parameters, Nmass):
     weight = float(Nmass - parameters.Nmin)/float(parameters.Nmax - parameters.Nmin)
     weight = min(1.0, weight)
     weight = max(0.0, weight)
-    
+
     return weight
 
-    
+
 
 def squared_distance(atom1, atom2):
 #    if atom1 in atom2.squared_distances:
@@ -518,9 +520,9 @@ def HydrogenBondEnergy(distance, dpka_max, cutoff, f_angle=1.0):
         value = 0.00
     else:
         value = 1.0-(distance-cutoff[0])/(cutoff[1]-cutoff[0])
-          
+
     dpKa  = dpka_max*value*f_angle
-    
+
     return abs(dpKa)
 
 
@@ -535,13 +537,13 @@ def AngleFactorX(atom1=None, atom2=None, atom3=None, center=None):
     dX_32 = atom2.x - atom3.x
     dY_32 = atom2.y - atom3.y
     dZ_32 = atom2.z - atom3.z
-    
+
     distance_23 = math.sqrt( dX_32*dX_32 + dY_32*dY_32 + dZ_32*dZ_32 )
-    
+
     dX_32 = dX_32/distance_23
     dY_32 = dY_32/distance_23
     dZ_32 = dZ_32/distance_23
-    
+
     if atom1 == None:
         dX_21 = center[0] - atom2.x
         dY_21 = center[1] - atom2.y
@@ -550,13 +552,13 @@ def AngleFactorX(atom1=None, atom2=None, atom3=None, center=None):
         dX_21 = atom1.x - atom2.x
         dY_21 = atom1.y - atom2.y
         dZ_21 = atom1.z - atom2.z
-        
+
     distance_12 = math.sqrt( dX_21*dX_21 + dY_21*dY_21 + dZ_21*dZ_21 )
-    
+
     dX_21 = dX_21/distance_12
     dY_21 = dY_21/distance_12
     dZ_21 = dZ_21/distance_12
-        
+
     f_angle = dX_21*dX_32 + dY_21*dY_32 + dZ_21*dZ_32
 
 
@@ -569,7 +571,7 @@ def hydrogen_bond_interaction(group1, group2, version):
     # find the smallest distance between interacting atoms
     atoms1 = group1.get_interaction_atoms(group2)
     atoms2 = group2.get_interaction_atoms(group1)
-    [closest_atom1, distance, closest_atom2] = Source.calculations.get_smallest_distance(atoms1, atoms2)
+    [closest_atom1, distance, closest_atom2] = propka.calculations.get_smallest_distance(atoms1, atoms2)
 
     if None in [closest_atom1, closest_atom2]:
         print('Warning: Side chain interaction failed for %s and %s'%(group1.label, group2.label))
@@ -577,20 +579,20 @@ def hydrogen_bond_interaction(group1, group2, version):
 
     # get the parameters
     [dpka_max, cutoff] = version.get_hydrogen_bond_parameters(closest_atom1,closest_atom2)
-    
+
     if dpka_max==None or None in cutoff:
         return None
-    
+
     # check that the closest atoms are close enough
     if distance >= cutoff[1]:
         return None
 
     # check that bond distance criteria is met
-    bond_distance_too_short = group1.atom.is_atom_within_bond_distance(group2.atom, 
+    bond_distance_too_short = group1.atom.is_atom_within_bond_distance(group2.atom,
                                                                        version.parameters.min_bond_distance_for_hydrogen_bonds,1)
     if bond_distance_too_short:
         return None
-        
+
     # set the angle factor
     #
     #  ---closest_atom1/2
@@ -602,18 +604,18 @@ def hydrogen_bond_interaction(group1, group2, version):
         if closest_atom2.element == 'H':
             heavy_atom = closest_atom2.bonded_atoms[0]
             hydrogen   = closest_atom2
-            distance, f_angle, nada = Source.calculations.AngleFactorX(closest_atom1, hydrogen, heavy_atom)   
+            distance, f_angle, nada = propka.calculations.AngleFactorX(closest_atom1, hydrogen, heavy_atom)
         else:
             # Either the structure is corrupt (no hydrogen), or the heavy atom is closer to
             # the titratable atom than the hydrogen. In either case we set the angle factor
             # to 0
             f_angle = 0.0
-   
+
     elif group1.type in version.parameters.angular_dependent_sidechain_interactions:
         if closest_atom1.element == 'H':
             heavy_atom = closest_atom1.bonded_atoms[0]
             hydrogen   = closest_atom1
-            distance, f_angle, nada = Source.calculations.AngleFactorX(closest_atom2, hydrogen, heavy_atom)  
+            distance, f_angle, nada = propka.calculations.AngleFactorX(closest_atom2, hydrogen, heavy_atom)
         else:
             # Either the structure is corrupt (no hydrogen), or the heavy atom is closer to
             # the titratable atom than the hydrogen. In either case we set the angle factor
@@ -621,7 +623,7 @@ def hydrogen_bond_interaction(group1, group2, version):
             f_angle = 0.0
 
     weight = version.calculatePairWeight(group1.Nmass, group2.Nmass)
-    
+
     exception, value = version.checkExceptions(group1, group2)
     #exception = False # circumventing exception
     if exception == True:
@@ -629,7 +631,7 @@ def hydrogen_bond_interaction(group1, group2, version):
             #print(" exception for %s %s %6.2lf" % (group1.label, group2.label, value))
     else:
         value = version.calculateSideChainEnergy(distance, dpka_max, cutoff, weight, f_angle)
-        
+
     # print('distance',distance)
     # print('dpka_max',dpka_max)
     # print('cutoff',cutoff)
@@ -654,7 +656,7 @@ def HydrogenBondEnergy(distance, dpka_max, cutoff, f_angle=1.0):
         value = 1.0-(distance-cutoff[0])/(cutoff[1]-cutoff[0])
 
     dpKa  = dpka_max*value*f_angle
-    
+
     return abs(dpKa)
 
 
@@ -678,15 +680,15 @@ def checkCoulombPair(parameters, group1, group2, distance):
     """
     Npair = group1.Nmass + group2.Nmass
     do_coulomb = True
-    
+
     # check if both groups are titratable (ions are taken care of in determinants::setIonDeterminants)
     if not (group1.titratable and group2.titratable):
         do_coulomb = False
-    
+
     # check if the distance is not too big
     if distance > parameters.coulomb_cutoff2:
         do_coulomb = False
-    
+
     # check that desolvation is ok
     if Npair < parameters.Nmin:
         do_coulomb = False
@@ -725,8 +727,8 @@ def BackBoneReorganization(parameters, conformation):
         dpKa = 0.00
         for BBC_group in BBC_groups:
             center = [titratable_group.x, titratable_group.y, titratable_group.z]
-            distance, f_angle, nada = AngleFactorX(atom2=BBC_group.get_interaction_atoms(titratable_group)[0], 
+            distance, f_angle, nada = AngleFactorX(atom2=BBC_group.get_interaction_atoms(titratable_group)[0],
-                                                   atom3=BBC_group.atom, 
+                                                   atom3=BBC_group.atom,
                                                    center=center)
             if distance <  6.0 and f_angle > 0.001:
                 value = 1.0-(distance-3.0)/(6.0-3.0)
@@ -767,7 +769,7 @@ def checkExceptions(version, group1, group2):
     else:
         # do nothing, no exception for this pair
         exception = False; value = None
-        
+
     return exception, value
 
 
@@ -780,7 +782,7 @@ def checkCooArgException(group_coo, group_arg, version):
     str = "xxx"
     exception = True
     value_tot = 0.00
-   
+
     #dpka_max = parameters.sidechain_interaction.get_value(group_coo.type,group_arg.type)
     #cutoff   = parameters.sidechain_cutoffs.get_value(group_coo.type,group_arg.type)
 
@@ -801,7 +803,7 @@ def checkCooArgException(group_coo, group_arg, version):
         if group_arg.type in version.parameters.angular_dependent_sidechain_interactions:
             atom3 = closest_arg_atom.bonded_atoms[0]
             distance, f_angle, nada = AngleFactorX(closest_coo_atom, closest_arg_atom, atom3)
-            
+
         value = HydrogenBondEnergy(distance, dpka_max, cutoff, f_angle)
         #print(iter, closest_coo_atom, closest_arg_atom,distance,value)
         value_tot += value
@@ -809,7 +811,7 @@ def checkCooArgException(group_coo, group_arg, version):
         atoms_coo.remove(closest_coo_atom)
         atoms_arg.remove(closest_arg_atom)
 
-    
+
     return exception, value_tot
 
 
@@ -820,7 +822,7 @@ def checkCooCooException(group1, group2, version):
     exception = True
     [closest_atom1, distance, closest_atom2] = get_smallest_distance(group1.get_interaction_atoms(group2),
                                                                      group2.get_interaction_atoms(group1))
-                                                                           
+
     #dpka_max = parameters.sidechain_interaction.get_value(group1.type,group2.type)
     #cutoff   = parameters.sidechain_cutoffs.get_value(group1.type,group2.type)
     [dpka_max, cutoff] = version.get_hydrogen_bond_parameters(closest_atom1,closest_atom2)
@@ -863,7 +865,7 @@ def checkCysHisException(group1, group2, version):
     if checkBuried(group1.Nmass, group2.Nmass):
         exception = True
 
-    return exception, version.parameters.CYS_HIS_exception 
+    return exception, version.parameters.CYS_HIS_exception
 
 
 def checkCysCysException(group1, group2, version):

propka/conformation_container.py

@@ -2,7 +2,10 @@
 # Container for molecular conformations
 #
 
-import Source.group, Source.determinants, Source.determinant, Source.ligand, Source.output, Source.coupled_groups, functools
+from __future__ import division
+from __future__ import print_function
+
+import propka.group, propka.determinants, propka.determinant, propka.ligand, propka.output, propka.coupled_groups, functools
 
 class Conformation_container:
     def __init__(self, name='', parameters=None, molecular_container=None):
@@ -13,10 +16,10 @@ class Conformation_container:
         self.groups = []
         self.chains = []
         self.current_iter_item = 0
-        
+
         self.marvin_pkas_calculated = False
         self.non_covalently_coupled_groups = False
-       
+
         return
 
 
@@ -28,21 +31,21 @@ class Conformation_container:
         for atom in self.get_non_hydrogen_atoms():
             # has this atom been checked for groups?
             if atom.groups_extracted == 0:
-                self.validate_group(Source.group.is_group(self.parameters, atom))
+                self.validate_group(propka.group.is_group(self.parameters, atom))
-            # if the atom has been checked in a another conformation, check if it has a 
+            # if the atom has been checked in a another conformation, check if it has a
-            # group that should be used in this conformation as well 
+            # group that should be used in this conformation as well
             elif atom.group:
                 self.validate_group(atom.group)
 
-        return 
+        return
 
     def additional_setup_when_reading_input_file(self):
-        
+
-        # if a group is coupled and we are reading a .propka_input file, 
+        # if a group is coupled and we are reading a .propka_input file,
         # some more configuration might be needed
 
         # print coupling map
-        map = Source.output.make_interaction_map('Covalent coupling map for %s'%self,
+        map = propka.output.make_interaction_map('Covalent coupling map for %s'%self,
                                                  self.get_covalently_coupled_groups(),
                                                  lambda g1,g2: g1 in g2.covalently_coupled_groups)
         print(map)
@@ -50,11 +53,11 @@ class Conformation_container:
         # check if we should set a common charge centre as well
         if self.parameters.common_charge_centre:
             self.set_common_charge_centres()
-                
+
         return
-    
+
     def set_common_charge_centres(self):
-        for system in self.get_coupled_systems(self.get_covalently_coupled_groups(), Source.group.Group.get_covalently_coupled_groups):
+        for system in self.get_coupled_systems(self.get_covalently_coupled_groups(), propka.group.Group.get_covalently_coupled_groups):
             # make a list of the charge centre coordinates
             all_coordinates = list(map(lambda g: [g.x, g.y, g.z], system))
             # find the common charge center
@@ -82,19 +85,19 @@ class Conformation_container:
                     continue
                 if cg.atom.sybyl_type == group.atom.sybyl_type:
                     group.couple_covalently(cg)
-                    
+
         # check if we should set a common charge centre as well
         if self.parameters.common_charge_centre:
             self.set_common_charge_centres()
 
         # print coupling map
-        map = Source.output.make_interaction_map('Covalent coupling map for %s'%self,
+        map = propka.output.make_interaction_map('Covalent coupling map for %s'%self,
                                                  #self.get_titratable_groups(),
                                                  self.get_covalently_coupled_groups(),
                                                  lambda g1,g2: g1 in g2.covalently_coupled_groups)
         print(map)
 
-            
+
         return
 
 
@@ -102,8 +105,8 @@ class Conformation_container:
         # check if non-covalent coupling has already been set up in an input file
         if len(list(filter(lambda g: len(g.non_covalently_coupled_groups)>0, self.get_titratable_groups())))>0:
             self.non_covalently_coupled_groups = True
-            
+
-        Source.coupled_groups.nccg.identify_non_covalently_coupled_groups(self,verbose=verbose)
+        propka.coupled_groups.nccg.identify_non_covalently_coupled_groups(self,verbose=verbose)
 
         # re-do the check
         if len(list(filter(lambda g: len(g.non_covalently_coupled_groups)>0, self.get_titratable_groups())))>0:
@@ -157,26 +160,26 @@ class Conformation_container:
 
     def calculate_pka(self, version, options):
         print('\nCalculating pKas for',self)
-                
+
         # calculate desolvation
         for group in self.get_titratable_groups()+self.get_ions():
             version.calculate_desolvation(group)
 
         # calculate backbone interactions
-        Source.determinants.setBackBoneDeterminants(self.get_titratable_groups(), self.get_backbone_groups(), version)
+        propka.determinants.setBackBoneDeterminants(self.get_titratable_groups(), self.get_backbone_groups(), version)
 
         # setting ion determinants
-        Source.determinants.setIonDeterminants(self, version)
+        propka.determinants.setIonDeterminants(self, version)
 
         # calculating the back-bone reorganization/desolvation term
         version.calculateBackBoneReorganization(self)
 
         # setting remaining non-iterative and iterative side-chain & Coulomb interaction determinants
-        Source.determinants.setDeterminants(self.get_sidechain_groups(), version=version, options=options)
+        propka.determinants.setDeterminants(self.get_sidechain_groups(), version=version, options=options)
 
         # calculating the total pKa values
         for group in self.groups: group.calculate_total_pka()
-           
+
         # take coupling effects into account
         penalised_labels = self.coupling_effects()
 
@@ -195,20 +198,20 @@ class Conformation_container:
 
     def coupling_effects(self):
         #
-        # Bases: The group with the highest pKa (the most stable one in the 
+        # Bases: The group with the highest pKa (the most stable one in the
-        #        charged form) will be the first one to adopt a proton as pH 
+        #        charged form) will be the first one to adopt a proton as pH
         #        is lowered and this group is allowed to titrate.
         #        The remaining groups are penalised
         #
         # Acids: The group with the highest pKa (the least stable one in the
-        #        charged form) will be the last group to loose the proton as 
+        #        charged form) will be the last group to loose the proton as
         #        pH is raised and will be penalised.
         #        The remaining groups are allowed to titrate.
         #
         penalised_labels = []
-        
+
-        for all_groups in self.get_coupled_systems(self.get_covalently_coupled_groups(), 
+        for all_groups in self.get_coupled_systems(self.get_covalently_coupled_groups(),
-                                                   Source.group.Group.get_covalently_coupled_groups):
+                                                   propka.group.Group.get_covalently_coupled_groups):
 
             # check if we should share determinants
             if self.parameters.shared_determinants:
@@ -227,14 +230,14 @@ class Conformation_container:
                 for a in all_groups:
                     if a == first_group:
                         continue # group with the highest pKa is allowed to titrate...
-                    a.coupled_titrating_group = first_group 
+                    a.coupled_titrating_group = first_group
                     penalised_labels.append(a.label) #... and the rest is penalised
 
         return penalised_labels
 
 
     def share_determinants(self, groups):
-        
+
         # make a list of the determinants to share
         types = ['sidechain','backbone','coulomb']
         for type in types:
@@ -250,7 +253,7 @@ class Conformation_container:
 
              # overwrite/add maximum value for each determinant
             for det_group in max_dets.keys():
-                new_determinant = Source.determinant.Determinant(det_group, max_dets[det_group])
+                new_determinant = propka.determinant.Determinant(det_group, max_dets[det_group])
                 for g in groups:
                     g.set_determinant(new_determinant,type)
 
@@ -260,7 +263,7 @@ class Conformation_container:
 
     def get_coupled_systems(self, groups, get_coupled_groups):
         """ This generator will yield one covalently coupled system at the time """
-        groups = set(groups)        
+        groups = set(groups)
         while len(groups)>0:
             # extract a system of coupled groups ...
             system = set()
@@ -272,7 +275,7 @@ class Conformation_container:
 
         return
 
-    
+
     def get_a_coupled_system_of_groups(self, new_group, coupled_groups, get_coupled_groups):
         coupled_groups.add(new_group)
         [self.get_a_coupled_system_of_groups(c, coupled_groups, get_coupled_groups) for c in get_coupled_groups(new_group) if c not in coupled_groups]
@@ -311,7 +314,7 @@ class Conformation_container:
         """ returns all sidechain groups needed for the pKa calculations """
         return [group for group in self.groups if ('BB' not in group.type\
                     and not group.atom.cysteine_bridge)]
-    
+
     def get_covalently_coupled_groups(self):
         return [g for g in self.groups if len(g.covalently_coupled_groups)>0]
 
@@ -331,7 +334,7 @@ class Conformation_container:
 
     def get_titratable_groups(self):
         return [group for group in self.groups if group.titratable]
-    
+
     def get_titratable_groups_and_cysteine_bridges(self):
         return [group for group in self.groups if group.titratable or group.residue_type == 'CYS']
 
@@ -367,19 +370,19 @@ class Conformation_container:
             atom.conformation_container = self
         if not atom.molecular_container:
             atom.molecular_container = self.molecular_container
-        
+
         # store chain id for bookkeeping
         if not atom.chainID in self.chains:
             self.chains.append(atom.chainID)
-            
+
-        return 
+        return
 
     def copy_atom(self, atom):
         new_atom  = atom.makeCopy()
         self.atoms.append(new_atom)
         new_atom.conformation_container = self
 
-        return 
+        return
 
     def get_non_hydrogen_atoms(self):
         return [atom for atom in self.atoms if atom.element!='H']
@@ -391,7 +394,7 @@ class Conformation_container:
             if not self.have_atom(atom):
                 self.copy_atom(atom)
         return
-    
+
     def have_atom(self, atom):
         res = False
         for a in self.atoms:
@@ -407,12 +410,12 @@ class Conformation_container:
                     return g
         return False
 
-    
+
     def set_ligand_atom_names(self):
         for atom in self.get_ligand_atoms():
-            Source.ligand.assign_sybyl_type(atom)
+            propka.ligand.assign_sybyl_type(atom)
         return
-            
+
 
 
     def __str__(self):
@@ -431,7 +434,7 @@ class Conformation_container:
         return
 
     def sort_atoms_key(self, atom):
-        key = ord(atom.chainID)*1e7 
+        key = ord(atom.chainID)*1e7
         key += atom.resNumb*1000
         if len(atom.name) > len(atom.element):
             key += ord(atom.name[len(atom.element)])

propka/coupled_groups.py

@@ -1,12 +1,15 @@
 
-import math, Source.output, Source.group, Source.lib, itertools
+from __future__ import division
+from __future__ import print_function
+
+import math, propka.output, propka.group, propka.lib, itertools
 
 
 class non_covalently_couple_groups:
     def __init__(self):
 
         self.parameters = None
-        
+
         # self.do_intrinsic = False
         # self.do_pair_wise = False
         self.do_prot_stat = True
@@ -18,10 +21,10 @@ class non_covalently_couple_groups:
     # Methods for finding coupled groups
     #
     def is_coupled_protonation_state_probability(self, group1, group2, energy_method, return_on_fail=True):
-    
+
         # check if the interaction energy is high enough
         interaction_energy = max(self.get_interaction(group1,group2), self.get_interaction(group2,group1))
- 
+
         if interaction_energy<=self.parameters.min_interaction_energy and return_on_fail:
             return {'coupling_factor':-1.0}
 
@@ -29,7 +32,7 @@ class non_covalently_couple_groups:
         for group in [group1, group2]:
             if not hasattr(group, 'intrinsic_pKa'):
                 group.calculate_intrinsic_pka()
-    
+
         use_pH = self.parameters.pH
         if self.parameters.pH == 'variable':
             use_pH = min(group1.pka_value, group2.pka_value)
@@ -61,7 +64,7 @@ class non_covalently_couple_groups:
         self.swap_interactions([group1], [group2])
         group1.calculate_total_pka()
         group2.calculate_total_pka()
-               
+
         # check difference in free energy
         if abs(default_energy - swapped_energy) > self.parameters.max_free_energy_diff and return_on_fail:
             return {'coupling_factor':-1.0}
@@ -80,21 +83,21 @@ class non_covalently_couple_groups:
             self.get_interaction_factor(interaction_energy)
 
         return {'coupling_factor':factor,
-                'default_energy':default_energy, 
+                'default_energy':default_energy,
-                'swapped_energy':swapped_energy, 
+                'swapped_energy':swapped_energy,
                 'interaction_energy':interaction_energy,
-                'swapped_pka1':swapped_pka1, 
+                'swapped_pka1':swapped_pka1,
                 'swapped_pka2':swapped_pka2,
                 'pka_shift1':pka_shift1,
                 'pka_shift2':pka_shift2,
                 'pH':use_pH}
 
 
-    
+
     #
     # Methods for calculating the coupling factor
-    # 
+    #
-    
+
     def get_pka_diff_factor(self, pka1, pka2):
         intrinsic_pka_diff = abs(pka1-pka2)
         res = 0.0
@@ -120,10 +123,10 @@ class non_covalently_couple_groups:
 
 
 
-    
+
     def identify_non_covalently_coupled_groups(self, conformation, verbose=True):
         """ Finds coupled residues in protein """
-    
+
         self.parameters = conformation.parameters
         if verbose:
             print('')
@@ -156,19 +159,19 @@ class non_covalently_couple_groups:
                         data = self.is_coupled_protonation_state_probability(g1, g2,conformation.calculate_folding_energy)
                         if data['coupling_factor'] >0.0:
                             g1.couple_non_covalently(g2)
-            
+
         if verbose:
             self.print_out_swaps(conformation)
 
         return
 
     def print_out_swaps(self, conformation, verbose=True):
-        map = Source.output.make_interaction_map('Non-covalent coupling map for %s'%conformation,
+        map = propka.output.make_interaction_map('Non-covalent coupling map for %s'%conformation,
                                                  conformation.get_non_covalently_coupled_groups(),
                                                  lambda g1,g2: g1 in g2.non_covalently_coupled_groups)
         print(map)
 
-        for system in conformation.get_coupled_systems(conformation.get_non_covalently_coupled_groups(),Source.group.Group.get_non_covalently_coupled_groups):
+        for system in conformation.get_coupled_systems(conformation.get_non_covalently_coupled_groups(),propka.group.Group.get_non_covalently_coupled_groups):
             self.print_system(conformation, list(system))
         return
 
@@ -187,7 +190,7 @@ class non_covalently_couple_groups:
         print(coup_info)
 
         # make list of possible combinations of swap to try out
-        combinations = Source.lib.generate_combinations(interactions)
+        combinations = propka.lib.generate_combinations(interactions)
 
         # Make possible swap combinations
         swap_info = ''
@@ -199,7 +202,7 @@ class non_covalently_couple_groups:
             for interaction in combination:
                 swap_info += ' %s %s\n'%(interaction[0].label, interaction[1].label)
 
-            # swap... 
+            # swap...
             for interaction in combination:
                 self.swap_interactions([interaction[0]],[interaction[1]])
 
@@ -213,13 +216,13 @@ class non_covalently_couple_groups:
         return
     #
     # Interaction and swapping methods
-    #                
+    #
 
     def get_interaction(self, group1, group2, include_side_chain_hbs = True):
         determinants = group1.determinants['coulomb']
         if include_side_chain_hbs:
             determinants = group1.determinants['sidechain'] + group1.determinants['coulomb']
-    
+
         interaction_energy = 0.0
         for det in determinants:
             if group2 == det.group:
@@ -234,7 +237,7 @@ class non_covalently_couple_groups:
         s = ' '+'-'*113+'\n'
         for group in system:
             s += self.tagged_print(' %-8s|'%tag,group.getDeterminantString(), all_labels)
-            
+
         return s+'\n'
 
     def swap_interactions(self, groups1, groups2, include_side_chain_hbs = True):
@@ -242,7 +245,7 @@ class non_covalently_couple_groups:
         for i in range(len(groups1)):
             group1 = groups1[i]
             group2 = groups2[i]
-            
+
             # swap the interactions!
             self.transfer_determinant(group1.determinants['coulomb'], group2.determinants['coulomb'], group1.label, group2.label)
             if include_side_chain_hbs:
@@ -266,7 +269,7 @@ class non_covalently_couple_groups:
         for det in determinants2:
             if det.label == label1:
                 from2to1.append(det)
-        
+
         # ...and transfer it!
         for det in from1to2:
             det.label = label1
@@ -280,7 +283,7 @@ class non_covalently_couple_groups:
 
         return
 
-    #  
+    #
     # Output methods
     #
 
@@ -295,7 +298,7 @@ class non_covalently_couple_groups:
     def make_data_to_string(self, data, group1, group2):
         s = """ %s and %s coupled (prot.state): %5.2f
  Energy levels:       %6.2f, %6.2f  (difference: %6.2f) at pH %6.2f
- Interaction energy:  %6.2f 
+ Interaction energy:  %6.2f
  Intrinsic pka's:     %6.2f, %6.2f  (difference: %6.2f)
  Swapped pKa's:       %6.2f, %6.2f  (difference: %6.2f, %6.2f)"""%(group1.label,
                                                                    group2.label,

propka/determinant.py

@@ -1,4 +1,7 @@
 
+from __future__ import division
+from __future__ import print_function
+
 class Determinant:
     """
         Determinant class - set up for later structurization

propka/determinants.py

@@ -1,8 +1,12 @@
+
+from __future__ import division
+from __future__ import print_function
+
 import math, time
 
-import Source.iterative, Source.lib, Source.vector_algebra
+import propka.iterative, propka.lib, propka.vector_algebra
-import Source.calculations
+import propka.calculations
-from Source.determinant import Determinant
+from propka.determinant import Determinant
 
 
 def setDeterminants(propka_groups, version=None, options=None):
@@ -19,19 +23,19 @@ def setDeterminants(propka_groups, version=None, options=None):
             # do not calculate interactions for coupled groups
             if group2 in group1.covalently_coupled_groups:
                 break
-            
+
-            distance = Source.calculations.distance(group1, group2)
+            distance = propka.calculations.distance(group1, group2)
-            
+
             if distance < version.parameters.coulomb_cutoff2:
                 interaction_type = version.parameters.interaction_matrix.get_value(group1.type,group2.type)
                 if interaction_type == 'I':
-                    Source.iterative.addtoDeterminantList(group1, group2, distance, iterative_interactions, version=version)
+                    propka.iterative.addtoDeterminantList(group1, group2, distance, iterative_interactions, version=version)
                 elif interaction_type == 'N':
                     addDeterminants(group1, group2, distance, version)
 
- 
+
     # --- Iterative section ---#
-    Source.iterative.addDeterminants(iterative_interactions, version, options=options)
+    propka.iterative.addDeterminants(iterative_interactions, version, options=options)
 
 
 def addDeterminants(group1, group2, distance, version):
@@ -52,11 +56,11 @@ def addSidechainDeterminants(group1, group2, version=None):
     adding side-chain determinants/perturbations
     Note, resNumb1 > resNumb2
     """
-    
+
     hbond_interaction = version.hydrogen_bond_interaction(group1, group2)
 
     if hbond_interaction:
-        
+
         if group1.charge == group2.charge:
             # acid pair or base pair
             if group1.model_pka < group2.model_pka:
@@ -78,7 +82,7 @@ def addCoulombDeterminants(group1, group2, distance, version):
     """
     adding NonIterative Coulomb determinants/perturbations
     """
-    
+
     coulomb_interaction = version.electrostatic_interaction(group1, group2, distance)
 
     if coulomb_interaction:
@@ -104,7 +108,7 @@ def addCoulombAcidPair(object1, object2, value):
     Adding the Coulomb interaction (an acid pair):
     the higher pKa is raised
     """
-    
+
     if object1.model_pka > object2.model_pka:
         newDeterminant = Determinant(object2, value)
         object1.determinants['coulomb'].append(newDeterminant)
@@ -151,7 +155,7 @@ def setIonDeterminants(conformation_container, version):
     """
     for titratable_group in conformation_container.get_titratable_groups():
         for ion_group in conformation_container.get_ions():
-            squared_distance = Source.calculations.squared_distance(titratable_group, ion_group)
+            squared_distance = propka.calculations.squared_distance(titratable_group, ion_group)
             if squared_distance < version.parameters.coulomb_cutoff2_squared:
                 weight = version.calculatePairWeight(titratable_group.Nmass, ion_group.Nmass)
                 # the pKa of both acids and bases are shifted up by negative ions (and vice versa)
@@ -162,16 +166,16 @@ def setIonDeterminants(conformation_container, version):
     return
 
 def setBackBoneDeterminants(titratable_groups, backbone_groups, version):
-    
+
     for titratable_group in titratable_groups:
         # find out which backbone groups this titratable is interacting with
         for backbone_group in backbone_groups:
             # find the interacting atoms
             backbone_interaction_atoms = backbone_group.get_interaction_atoms(titratable_group)
             titratable_group_interaction_atoms = titratable_group.interaction_atoms_for_acids
-            
+
             # find the smallest distance
-            [backbone_atom, distance, titratable_atom] = Source.calculations.get_smallest_distance(backbone_interaction_atoms,
+            [backbone_atom, distance, titratable_atom] = propka.calculations.get_smallest_distance(backbone_interaction_atoms,
                                                                                                    titratable_group_interaction_atoms)
             # get the parameters
             parameters = version.get_backbone_hydrogen_bond_parameters(backbone_atom, titratable_atom)
@@ -186,7 +190,7 @@ def setBackBoneDeterminants(titratable_groups, backbone_groups, version):
                 # for BBC groups, the hydrogen is on the titratable group
                 #
                 #        Titra.
-                #       / 
+                #       /
                 #      H
                 #     .
                 #    O
@@ -197,7 +201,7 @@ def setBackBoneDeterminants(titratable_groups, backbone_groups, version):
                         if titratable_atom.element == 'H':
                             heavy_atom    = titratable_atom.bonded_atoms[0]
                             hydrogen_atom = titratable_atom
-                            [d1, f_angle, d2] = Source.calculations.AngleFactorX(atom1=heavy_atom,
+                            [d1, f_angle, d2] = propka.calculations.AngleFactorX(atom1=heavy_atom,
                                                                                  atom2=hydrogen_atom,
                                                                                  atom3=backbone_atom)
                         else:
@@ -218,7 +222,7 @@ def setBackBoneDeterminants(titratable_groups, backbone_groups, version):
                     if backbone_atom.element == 'H':
                         backbone_N = backbone_atom.bonded_atoms[0]
                         backbone_H = backbone_atom
-                        [d1, f_angle, d2] = Source.calculations.AngleFactorX(atom1=titratable_atom,
+                        [d1, f_angle, d2] = propka.calculations.AngleFactorX(atom1=titratable_atom,
                                                                              atom2=backbone_H,
                                                                              atom3=backbone_N)
                     else:
@@ -226,13 +230,13 @@ def setBackBoneDeterminants(titratable_groups, backbone_groups, version):
                         # the titratable atom than the hydrogen. In either case we set the angle factor
                         # to 0
                         f_angle = 0.0
-                    
+
-            
+
                 if f_angle > 0.001:
-                    value = titratable_group.charge * Source.calculations.HydrogenBondEnergy(distance, dpKa_max, [cutoff1,cutoff2], f_angle)
+                    value = titratable_group.charge * propka.calculations.HydrogenBondEnergy(distance, dpKa_max, [cutoff1,cutoff2], f_angle)
-                    
+
                     newDeterminant = Determinant(backbone_group, value)
                     titratable_group.determinants['backbone'].append(newDeterminant)
 
-        
+
     return

propka/group.py

@@ -2,9 +2,12 @@
 # Class for storing groups important for propka calculations
 #
 
-import Source.ligand, Source.determinant, Source.ligand_pka_values, math, Source.protonate
+from __future__ import division
+from __future__ import print_function
 
-my_protonator = Source.protonate.Protonate(verbose=False)
+import propka.ligand, propka.determinant, propka.ligand_pka_values, math, propka.protonate
+
+my_protonator = propka.protonate.Protonate(verbose=False)
 
 expected_atoms_acid_interactions = {
     'COO':{'O':2},
@@ -65,7 +68,7 @@ expected_atoms_base_interactions = {
     'SH': {'S':1},
     'CG': {'N':3},
     'C2N':{'N':2},
-    'OCO':{'O':2},  
+    'OCO':{'O':2},
     'N30':{'N':1},
     'N31':{'N':1},
     'N32':{'N':1},
@@ -86,7 +89,7 @@ class Group:
         self.determinants = {'sidechain':[],'backbone':[],'coulomb':[]}
         self.pka_value = 0.0
         self.model_pka = 0.0
-        
+
         self.Emass = 0.0
         self.Nmass = 0.0
         self.Elocl = 0.0
@@ -99,18 +102,18 @@ class Group:
         self.interaction_atoms_for_acids = []
         self.interaction_atoms_for_bases = []
         self.model_pka_set = False
-        
+
         # information on covalent and non-covalent coupling
         self.non_covalently_coupled_groups = []
         self.covalently_coupled_groups = []
         self.coupled_titrating_group = None
         self.common_charge_centre = False
-        
+
 
         self.residue_type = self.atom.resName
         if self.atom.terminal:
             self.residue_type = self.atom.terminal
-       
+
         if self.atom.type=='atom':
             self.label = '%-3s%4d%2s'%(self.residue_type, atom.resNumb, atom.chainID)
         elif self.atom.resName in ['DA ','DC ','DG ','DT ']:
@@ -119,7 +122,7 @@ class Group:
                                             atom.name.replace('\'','')[-1],
                                             atom.resNumb,
                                             atom.chainID)
-            
+
 #            self.label = '%1s%1s%1s%4d%2s'%(self.residue_type[1], atom.element,atom.name[-1], atom.resNumb, atom.chainID)
         else:
             self.label = '%-3s%4s%2s'%(self.residue_type, atom.name, atom.chainID)
@@ -133,28 +136,28 @@ class Group:
 
     #
     # Coupling-related methods
-    #   
+    #
     def couple_covalently(self, other):
         """ Couple this group with another group """
         # do the coupling
         if not other in self.covalently_coupled_groups:
             self.covalently_coupled_groups.append(other)
-        
+
         if not self in other.covalently_coupled_groups:
             other.covalently_coupled_groups.append(self)
 
-        return 
+        return
 
     def couple_non_covalently(self, other):
         """ Couple this group with another group """
         # do the coupling
         if not other in self.non_covalently_coupled_groups:
             self.non_covalently_coupled_groups.append(other)
-        
+
         if not self in other.non_covalently_coupled_groups:
             other.non_covalently_coupled_groups.append(self)
 
-        return 
+        return
 
     def get_covalently_coupled_groups(self): return self.covalently_coupled_groups
     def get_non_covalently_coupled_groups(self): return self.non_covalently_coupled_groups
@@ -168,12 +171,12 @@ class Group:
                 continue
 
             for type in ['sidechain','backbone','coulomb']:
-                for g in other.determinants[type]: self.share_determinant(g,type) 
+                for g in other.determinants[type]: self.share_determinant(g,type)
-                  
+
         # recalculate pka values
         self.calculate_total_pka()
         other.calculate_total_pka()
-            
+
         return
 
 
@@ -190,19 +193,19 @@ class Group:
 
         # otherwise we just add the determinant to our list
         if not added:
-            self.determinants[type].append(Source.determinant.Determinant(new_determinant.group, 
+            self.determinants[type].append(propka.determinant.Determinant(new_determinant.group,
                                                                           new_determinant.value))
 
         return
 
     def make_covalently_coupled_line(self):
-        
+
         # first check if there are any coupled groups at all
         if len(self.covalently_coupled_groups) == 0:
             return ''
 
         line = 'CCOUPL%5d'%self.atom.numb
-        
+
         # extract and sort numbers of coupled groups
         coupled = []
         for g in self.covalently_coupled_groups:
@@ -213,7 +216,7 @@ class Group:
         for b in coupled:
             line += '%5d'%b
         line += '\n'
-       
+
         return line
 
     def make_non_covalently_coupled_line(self):
@@ -222,7 +225,7 @@ class Group:
             return ''
 
         line = 'NCOUPL%5d'%self.atom.numb
-        
+
         # extract and sort numbers of coupled groups
         coupled = []
         for g in self.non_covalently_coupled_groups:
@@ -233,7 +236,7 @@ class Group:
         for b in coupled:
             line += '%5d'%b
         line += '\n'
-    
+
 
         return line
 
@@ -283,7 +286,7 @@ class Group:
                 return
 
         # otherwise we just add the determinant to our list
-        self.determinants[type].append(Source.determinant.Determinant(new_determinant.group, 
+        self.determinants[type].append(propka.determinant.Determinant(new_determinant.group,
                                                                       new_determinant.value))
 
         return
@@ -298,7 +301,7 @@ class Group:
                 return
 
         # otherwise we just add the determinant to our list
-        self.determinants[type].append(Source.determinant.Determinant(new_determinant.group, 
+        self.determinants[type].append(propka.determinant.Determinant(new_determinant.group,
                                                                       new_determinant.value))
 
         return
@@ -344,12 +347,12 @@ class Group:
             self.charge = self.parameters.charge[self.type]
         if self.residue_type in self.parameters.ions.keys():
             self.charge = self.parameters.ions[self.residue_type]
-            
+
         #print('ION setup',self,self.residue_type, self.charge)
 
         # find the center and the interaction atoms
         self.setup_atoms()
-     
+
         # set the model pka value
         self.titratable = False
         if self.residue_type in self.parameters.model_pkas.keys():
@@ -365,7 +368,7 @@ class Group:
         if self.model_pka_set and not self.atom.cysteine_bridge:
             self.titratable = True
 
-    
+
         return
 
     def setup_atoms(self):
@@ -378,7 +381,7 @@ class Group:
 
     def set_interaction_atoms(self, interaction_atoms_for_acids, interaction_atoms_for_bases):
         [a.set_group_type(self.type) for a in interaction_atoms_for_acids+interaction_atoms_for_bases]
-                
+
         self.interaction_atoms_for_acids = interaction_atoms_for_acids
         self.interaction_atoms_for_bases = interaction_atoms_for_bases
 
@@ -400,11 +403,11 @@ class Group:
             print('         Expected %d interaction atoms for acids, found:'%Na)
             for i in range(len(self.interaction_atoms_for_acids)):
                  print('             %s'%self.interaction_atoms_for_acids[i])
-                
+
             print('         Expected %d interaction atoms for bases, found:'%Nb)
             for i in range(len(self.interaction_atoms_for_bases)):
                  print('             %s'%self.interaction_atoms_for_bases[i])
-                
+
 
                     #return
 
@@ -437,7 +440,7 @@ class Group:
         number_of_sidechain = len(self.determinants['sidechain'])
         number_of_backbone  = len(self.determinants['backbone'])
         number_of_coulomb   = len(self.determinants['coulomb'])
-        
+
         number_of_lines     = max(1, number_of_sidechain, number_of_backbone, number_of_coulomb)
         str  = ""
         for line_number in range(number_of_lines):
@@ -458,11 +461,11 @@ class Group:
                 str += " %6.2lf %4d" % (self.Elocl, self.Nlocl)
             else:
                 str += "%40s" % (" ")
-        
+
             # add the determinants
             for type in ['sidechain','backbone','coulomb']:
                 str += self.get_determinant_for_string(type,line_number)
-            
+
             # adding end-of-line
             str += "\n"
 
@@ -484,7 +487,7 @@ class Group:
         if self.atom.cysteine_bridge:
             self.pka_value = 99.99
             return
-        
+
 
         self.pka_value = self.model_pka + self.Emass + self.Elocl
 
@@ -494,12 +497,12 @@ class Group:
 
         return
 
-        
+
 
 
     def calculate_intrinsic_pka(self):
         """
-        Calculates the intrinsic pKa values from the desolvation determinants, back-bone hydrogen bonds, 
+        Calculates the intrinsic pKa values from the desolvation determinants, back-bone hydrogen bonds,
         and side-chain hydrogen bond to non-titratable residues
         """
         back_bone  = 0.0
@@ -528,18 +531,18 @@ class Group:
         ligand_type = ''
         if self.atom.type == 'hetatm':
             ligand_type = self.type
-        
+
         penalty = ''
         if self.coupled_titrating_group:
             penalty = ' NB: Discarded due to coupling with %s'%self.coupled_titrating_group.label
 
-        str = "   %9s %8.2lf %10.2lf %18s   %s\n" % (self.label, 
+        str = "   %9s %8.2lf %10.2lf %18s   %s\n" % (self.label,
-                                                     self.pka_value, 
+                                                     self.pka_value,
-                                                     self.model_pka,ligand_type, 
+                                                     self.model_pka,ligand_type,
                                                      penalty)
 
         return str
-     
+
     def __str__(self):
         return 'Group (%s) for %s'%(self.type,self.atom)
 
@@ -559,7 +562,7 @@ class Group:
             reference = parameters.reference
 
         # If not titratable, the contribution is zero
-        
+
         if not self.titratable:
             return 0.00
 
@@ -571,25 +574,25 @@ class Group:
                 if determinant.value > 0.00:
                     pka_prime -= determinant.value
             ddG_neutral = -1.36*(pka_prime - self.model_pka)
-                    
+
         # calculating the ddG(low-pH --> pH) contribution
         # folded
         x =  pH - self.pka_value
         y = 10**x
         Q_pro = math.log10(1+y)
-        
+
         # unfolded
         x =  pH - self.model_pka
         y = 10**x
         Q_mod = math.log10(1+y)
-        
+
         ddG_low = -1.36*(Q_pro - Q_mod)
         ddG = ddG_neutral + ddG_low
 
         return ddG
 
     def calculate_charge(self, parmaeters, pH=7.0, state='folded'):
-        
+
         if state == "unfolded":
             x =  self.charge * (self.model_pka - pH)
         else:
@@ -597,7 +600,7 @@ class Group:
 
         y = 10**x
         charge = self.charge*(y/(1.0+y))
-        
+
         return charge
 
 
@@ -610,7 +613,7 @@ class COO_group(Group):
     def setup_atoms(self):
         # Identify the two caroxyl oxygen atoms
         the_oxygens = self.atom.get_bonded_elements('O')
- 
+
         # set the center using the two oxygen carboxyl atoms
         self.set_center(the_oxygens)
         self.set_interaction_atoms(the_oxygens, the_oxygens)
@@ -624,7 +627,7 @@ class HIS_group(Group):
 
     def setup_atoms(self):
         # Find the atoms in the histidine ring
-        ring_atoms = Source.ligand.is_ring_member(self.atom)
+        ring_atoms = propka.ligand.is_ring_member(self.atom)
         if len(ring_atoms) != 5:
             print('Warning: His group does not seem to contain a ring',self)
 
@@ -634,11 +637,11 @@ class HIS_group(Group):
 
         # set the center using the ring atoms
         self.set_center(ring_atoms)
-        
+
         # find the hydrogens on the ring-nitrogens
         hydrogens = []
         nitrogens = [ra for ra in ring_atoms if ra.element == 'N']
-       
+
         for nitrogen in nitrogens:
             hydrogens.extend(nitrogen.get_bonded_elements('H'))
 
@@ -673,12 +676,12 @@ class ARG_group(Group):
     def setup_atoms(self):
         # set the center at the position of the main atom
         self.set_center([self.atom])
-        
+
         # find all the hydrogens on the nitrogen atoms
         nitrogens = self.atom.get_bonded_elements('N')
         for n in nitrogens:
             my_protonator.protonate_atom(n)
-        
+
         hydrogens = []
         for nitrogen in nitrogens:
             hydrogens.extend(nitrogen.get_bonded_elements('H'))
@@ -705,19 +708,19 @@ class AMD_group(Group):
     def setup_atoms(self):
         # Identify the oxygen and nitrogen amide atoms
         the_oxygen = self.atom.get_bonded_elements('O')
-        the_nitrogen = self.atom.get_bonded_elements('N') 
+        the_nitrogen = self.atom.get_bonded_elements('N')
-        
+
-        # add protons to the nitrogen 
+        # add protons to the nitrogen
         my_protonator.protonate_atom(the_nitrogen[0])
         the_hydrogens = the_nitrogen[0].get_bonded_elements('H')
 
         # set the center using the oxygen and nitrogen amide atoms
         self.set_center(the_oxygen+the_nitrogen)
-        
+
         self.set_interaction_atoms(the_nitrogen+the_hydrogens,the_oxygen)
 
         return
-    
+
 
 class TRP_group(Group):
     def __init__(self, atom):
@@ -727,7 +730,7 @@ class TRP_group(Group):
     def setup_atoms(self):
         # set the center at the position of the main atom
         self.set_center([self.atom])
-                 
+
         # find the hydrogen on the nitrogen atom
         my_protonator.protonate_atom(self.atom)
         the_hydrogen = self.atom.get_bonded_elements('H')
@@ -751,12 +754,12 @@ class Cterm_group(Group):
         the_carbon = self.atom.get_bonded_elements('C')
         the_other_oxygen = the_carbon[0].get_bonded_elements('O')
         the_other_oxygen.remove(self.atom)
-        
+
         # set the center and interaction atoms
         the_oxygens = [self.atom]+ the_other_oxygen
         self.set_center(the_oxygens)
         self.set_interaction_atoms(the_oxygens, the_oxygens)
-                
+
 
         return
 
@@ -802,7 +805,7 @@ class NAR_group(Group):
         self.residue_type = 'NAR'
         print('Found NAR group:',atom)
         return
-      
+
 
     def setup_atoms(self):
         # Identify the hydrogen
@@ -825,7 +828,7 @@ class NAM_group(Group):
         self.residue_type = 'NAM'
         print('Found NAM group:',atom)
         return
-        
+
 
     def setup_atoms(self):
         # Identify the hydrogen
@@ -936,7 +939,7 @@ class CG_group(Group):
 
         # set the center using the nitrogen
         self.set_center([self.atom])
-        
+
         the_hydrogens = []
         for n in the_nitrogens:
             my_protonator.protonate_atom(n)
@@ -1065,7 +1068,7 @@ class NP1_group(Group):
         print('Found NP1 group:',atom)
         return
 
- 
+
     def setup_atoms(self):
         # Identify the nitrogens
         my_protonator.protonate_atom(self.atom)
@@ -1123,16 +1126,16 @@ class titratable_ligand_group(Group):
             self.model_pka = atom.marvin_pka
             print('Titratable ligand group    ',atom, self.model_pka, self.charge)
         self.model_pka_set = True
-        
+
         return
 
 
 def is_group(parameters, atom):
     atom.groups_extracted = 1
-    
+
     # check if this atom belongs to a protein group
     protein_group = is_protein_group(parameters, atom)
-    if protein_group: return protein_group 
+    if protein_group: return protein_group
 
     # check if this atom belongs to a ion group
     ion_group = is_ion_group(parameters, atom)
@@ -1168,7 +1171,7 @@ def is_protein_group(parameters,atom):
     ### Backbone
     if atom.type == 'atom' and atom.name == 'N':
         # ignore proline backbone nitrogens
-        if atom.resName != 'PRO': 
+        if atom.resName != 'PRO':
             return BBN_group(atom)
     if atom.type == 'atom' and atom.name == 'C':
         # ignore C- carboxyl
@@ -1199,10 +1202,10 @@ def is_ligand_group_by_groups(parameters, atom):
     if atom.sybyl_type == 'N.ar':
         if len(atom.get_bonded_heavy_atoms())==2:
             return NAR_group(atom)
-    
+
     if atom.sybyl_type == 'N.am':
         return NAM_group(atom)
-    
+
     if atom.sybyl_type in ['N.3', 'N.4']:
         heavy_bonded = atom.get_bonded_heavy_atoms()
         if len(heavy_bonded) == 0:
@@ -1218,14 +1221,14 @@ def is_ligand_group_by_groups(parameters, atom):
         return N1_group(atom)
 
     if atom.sybyl_type == 'N.pl3':
-        # make sure that this atom is not part of a guadinium or amidinium group 
+        # make sure that this atom is not part of a guadinium or amidinium group
         bonded_carbons = atom.get_bonded_elements('C')
         if len(bonded_carbons) == 1:
             bonded_nitrogens = bonded_carbons[0].get_bonded_elements('N')
             if len(bonded_nitrogens) == 1:
                 return NP1_group(atom)
 
-        
+
     if atom.sybyl_type == 'C.2':
         # Guadinium and amidinium groups
         bonded_nitrogens = atom.get_bonded_elements('N')
@@ -1263,7 +1266,7 @@ def is_ligand_group_by_groups(parameters, atom):
 
     if atom.sybyl_type == 'O.2':
         return O2_group(atom)
-    
+
 
     if atom.sybyl_type == 'S.3':
         # sulphide group
@@ -1272,8 +1275,8 @@ def is_ligand_group_by_groups(parameters, atom):
             # other SP3 oxygen
         #else:
         #    return S3_group(atom)
-        
+
-    
+
     return None
 
 
@@ -1284,11 +1287,11 @@ def is_ligand_group_by_marvin_pkas(parameters, atom):
     # calculate Marvin ligand pkas for this conformation container
     # if not already done
     if not atom.conformation_container.marvin_pkas_calculated:
-        lpka = Source.ligand_pka_values.ligand_pka_values(parameters)
+        lpka = propka.ligand_pka_values.ligand_pka_values(parameters)
         lpka.get_marvin_pkas_for_molecular_container(atom.molecular_container,
                                                      min_pH=parameters.min_ligand_model_pka,
                                                      max_pH=parameters.max_ligand_model_pka)
-    
+
 
     if atom.marvin_pka:
         return titratable_ligand_group(atom)

propka/iterative.py

@@ -1,7 +1,12 @@
+
+from __future__ import division
+from __future__ import print_function
+
 import math, time
-import Source.lib as lib
+
-from Source.determinant import Determinant
+import propka.lib as lib
-import Source.calculations
+from propka.determinant import Determinant
+import propka.calculations
 
 # Some library functions for the interative pKa determinants
 
@@ -32,11 +37,11 @@ def addtoDeterminantList(group1, group2, distance, iterative_interactions, versi
 
 
 def addIterativeAcidPair(object1, object2, interaction):
-    """ 
+    """
     Adding the Coulomb 'iterative' interaction (an acid pair):
     the higher pKa is raised  with QQ+HB
     the lower  pKa is lowered with HB
-    """ 
+    """
     values       = interaction[1]
     annihilation = interaction[2]
     hbond_value   = values[0]
@@ -69,10 +74,10 @@ def addIterativeAcidPair(object1, object2, interaction):
 
 
 def addIterativeBasePair(object1, object2, interaction):
-    """ 
+    """
     Adding the Coulomb 'iterative' interaction (a base pair):
     the lower pKa is lowered
-    """ 
+    """
     values       = interaction[1]
     annihilation = interaction[2]
     hbond_value   = values[0]
@@ -106,10 +111,10 @@ def addIterativeBasePair(object1, object2, interaction):
 
 
 def addIterativeIonPair(object1, object2, interaction, version):
-    """ 
+    """
     Adding the Coulomb 'iterative' interaction (an acid-base pair):
     the pKa of the acid is lowered & the pKa of the base is raised
-    """ 
+    """
     values       = interaction[1]
     annihilation = interaction[2]
     hbond_value   = values[0]
@@ -132,7 +137,7 @@ def addIterativeIonPair(object1, object2, interaction, version):
 
     annihilation[0] = 0.00
     annihilation[1] = 0.00
-    
+
     if add_term == True:
 
       # Coulomb
@@ -161,9 +166,9 @@ def addIterativeIonPair(object1, object2, interaction, version):
 
 
 def addDeterminants(iterative_interactions, version, options=None):
-    """ 
+    """
     The iterative pKa scheme. Later it is all added in 'calculateTotalPKA'
-    """ 
+    """
     # --- setup ---
     iteratives = []
     done_group = []
@@ -270,7 +275,7 @@ def addDeterminants(iterative_interactions, version, options=None):
                     g = interaction[0]
                     newDeterminant = Determinant(g, value)
                     itres.group.determinants[type].append(newDeterminant)
-    
+
 
 
 def findIterative(pair, iteratives):

propka/lib.py

@@ -1,23 +1,27 @@
-#!/usr/bin/python
+from __future__ import division
+from __future__ import print_function
 
 import string, sys, copy, math, os
 
+import pkg_resources
 
 #
 # file I/O
 #
 def open_file_for_reading(filename):
-    if not os.path.isfile(filename):
+    try:
+        f = open(filename,'r')
+    except:
         raise Exception('Cannot find file %s' %filename)
-    
+    return f
-    return open(filename,'r')
 
 def open_file_for_writing(filename):
-    res = open(filename,'w')
+    try:
-    if not res:
+        res = open(filename,'w')
+    except:
         raise Exception('Could not open %s'%filename)
     return res
-    
+
 #
 # bookkeeping etc.
 #
@@ -43,7 +47,7 @@ def make_molecule(atom, atoms):
         if ba in atoms:
             atoms.remove(ba)
             res_atoms.extend(make_molecule(ba, atoms))
-         
+
     return res_atoms
 
 
@@ -61,7 +65,7 @@ def generate_combinations(interactions):
     res.remove([])
 
     return res
-    
+
 
 def make_combination(combis, interaction):
     res = []
@@ -85,37 +89,37 @@ def loadOptions():
     parser = OptionParser(usage)
 
     # loading the parser
-    parser.add_option("-f", "--file", action="append", dest="filenames", 
+    parser.add_option("-f", "--file", action="append", dest="filenames",
            help="read data from <filename>, i.e. <filename> is added to arguments")
-    parser.add_option("-r", "--reference", dest="reference", default="neutral", 
+    parser.add_option("-r", "--reference", dest="reference", default="neutral",
            help="setting which reference to use for stability calculations [neutral/low-pH]")
-    parser.add_option("-c", "--chain", action="append", dest="chains", 
+    parser.add_option("-c", "--chain", action="append", dest="chains",
            help="creating the protein with only a specified chain, note, chains without ID are labeled 'A' [all]")
-    parser.add_option("-t", "--thermophile", action="append", dest="thermophiles", 
+    parser.add_option("-t", "--thermophile", action="append", dest="thermophiles",
            help="defining a thermophile filename; usually used in 'alignment-mutations'")
-    parser.add_option("-a", "--alignment", action="append", dest="alignment", 
+    parser.add_option("-a", "--alignment", action="append", dest="alignment",
            help="alignment file connecting <filename> and <thermophile> [<thermophile>.pir]")
-    parser.add_option("-m", "--mutation", action="append", dest="mutations", 
+    parser.add_option("-m", "--mutation", action="append", dest="mutations",
            help="specifying mutation labels which is used to modify <filename> according to, e.g. N25R/N181D")
-    parser.add_option("-v", "--version", dest="version_label", default="Jan15", 
+    parser.add_option("-v", "--version", dest="version_label", default="Jan15",
            help="specifying the sub-version of propka [Jan15/Dec19]")
-    parser.add_option("-p", "--parameters",dest="parameters", default="propka.cfg",
+    parser.add_option("-p", "--parameters",dest="parameters", default=pkg_resources.resource_filename(__name__, "propka.cfg"),
-                      help="set the parameter file")
+                      help="set the parameter file [%default]")
-    parser.add_option("-z", "--verbose", dest="verbose", action="store_true", default=True, 
+    parser.add_option("-z", "--verbose", dest="verbose", action="store_true", default=True,
            help="sleep during calculations")
     parser.add_option("-q", "--quiet", dest="verbose", action="store_false",
            help="sleep during calculations")
-    parser.add_option("-s", "--silent",  dest="verbose", action="store_false", 
+    parser.add_option("-s", "--silent",  dest="verbose", action="store_false",
            help="not activated yet")
-    parser.add_option("--verbosity",  dest="verbosity", action="store_const", 
+    parser.add_option("--verbosity",  dest="verbosity", action="store_const",
            help="level of printout - not activated yet")
-    parser.add_option("-o", "--pH", dest="pH", type="float", default=7.0, 
+    parser.add_option("-o", "--pH", dest="pH", type="float", default=7.0,
            help="setting pH-value used in e.g. stability calculations [7.0]")
     parser.add_option("-w", "--window", dest="window", nargs=3, type="float", default=(0.0, 14.0, 1.0),
            help="setting the pH-window to show e.g. stability profiles [0.0, 14.0, 1.0]")
     parser.add_option("-g", "--grid",   dest="grid",   nargs=3, type="float", default=(0.0, 14.0, 0.1),
            help="setting the pH-grid to calculate e.g. stability related properties [0.0, 14.0, 0.1]")
-    parser.add_option("--mutator", dest="mutator", 
+    parser.add_option("--mutator", dest="mutator",
            help="setting approach for mutating <filename> [alignment/scwrl/jackal]")
     parser.add_option("--mutator-option", dest="mutator_options", action="append",
            help="setting property for mutator [e.g. type=\"side-chain\"]")
@@ -156,7 +160,7 @@ def makeTidyAtomLabel(name,element):
     """
     Returns a 'tidier' atom label for printing the new pdbfile
     """
-    
+
     if len(name)>4:# if longer than 4, just truncate the name
         label=name[0:4]
     elif len(name)==4:# if lenght is 4, otherwise use the name as it is

propka/ligand.py

@@ -1,7 +1,12 @@
 #!/usr/bin/python
 
-import sys, Source.calculations
+from __future__ import division
-from Source.vector_algebra import *
+from __future__ import print_function
+
+import sys
+
+import propka.calculations
+from propka.vector_algebra import *
 
 
 all_sybyl_types = [
@@ -157,7 +162,7 @@ def assign_sybyl_type(atom):
         bonded_elements[i]=atom.bonded_atoms[i].element
 
 
-    
+
     # Aromatic carbon/nitrogen
     if aromatic:
         for ra in ring_atoms:
@@ -187,7 +192,7 @@ def assign_sybyl_type(atom):
             nitrogens = atom.get_bonded_elements('N')
             oxygens = atom.get_bonded_elements('O')
             if len(nitrogens)==1 and len(oxygens)==1:
-                C = atom 
+                C = atom
                 N = nitrogens[0]
                 O = oxygens[0]
 
@@ -199,9 +204,9 @@ def assign_sybyl_type(atom):
                 set_type(C,'C.2')
                 set_type(O,'O.2')
                 return
-                
 
-    if atom.element=='C':                
+
+    if atom.element=='C':
         # check for carboxyl
         if len(atom.bonded_atoms)==3 and list(bonded_elements.values()).count('O')==2:
             i1 = list(bonded_elements.values()).index('O')
@@ -213,11 +218,11 @@ def assign_sybyl_type(atom):
                 return
 
 
-                
+
         # sp carbon
         if len(atom.bonded_atoms)<=2:
             for b in atom.bonded_atoms:
-                if Source.calculations.squared_distance(atom, b) < max_C_triple_bond_squared:
+                if propka.calculations.squared_distance(atom, b) < max_C_triple_bond_squared:
                     set_type(atom,'C.1')
                     set_type(b,b.element+'.1')
             if atom.sybyl_assigned:
@@ -232,7 +237,7 @@ def assign_sybyl_type(atom):
                     if len(b.bonded_atoms)<3 or is_planar(b):
                         set_type(b,'N.pl3')
             return
-        
+
         # sp3 carbon
         set_type(atom, 'C.3')
         return
@@ -244,11 +249,11 @@ def assign_sybyl_type(atom):
             if is_planar(atom.bonded_atoms[0]):
                 set_type(atom,'N.pl3')
                 return
-                    
+
         if planar:
             set_type(atom,'N.pl3')
             return
-                
+
         set_type(atom,'N.3')
         return
 
@@ -262,7 +267,7 @@ def assign_sybyl_type(atom):
                 the_carbon = atom.bonded_atoms[0]
                 if len(the_carbon.bonded_atoms)==3 and the_carbon.count_bonded_elements('O')==2:
                     [O1,O2] = the_carbon.get_bonded_elements('O')
-                    
+
                     if len(O1.bonded_atoms)==1 and len(O2.bonded_atoms)==1:
                         set_type(O1,'O.co2-')
                         set_type(O2,'O.co2')
@@ -270,7 +275,7 @@ def assign_sybyl_type(atom):
                         return
 
             # check for X=O
-            if Source.calculations.squared_distance(atom, atom.bonded_atoms[0]) < max_C_double_bond_squared:
+            if propka.calculations.squared_distance(atom, atom.bonded_atoms[0]) < max_C_double_bond_squared:
                 set_type(atom,'O.2')
                 if atom.bonded_atoms[0].element=='C':
                     set_type(atom.bonded_atoms[0],'C.2')
@@ -287,7 +292,7 @@ def assign_sybyl_type(atom):
             set_type(atom.bonded_atoms[i2],'O.2')
             set_type(atom,'S.o2')
             return
-                
+
         # check for SO4
         if list(bonded_elements.values()).count('O')==4:
             no_O2 = 0
@@ -297,7 +302,7 @@ def assign_sybyl_type(atom):
                     no_O2+=1
                 else:
                     set_type(atom.bonded_atoms[i],'O.3')
-                    
+
         set_type(atom,'S.3')
 
 
@@ -315,8 +320,8 @@ def assign_sybyl_type(atom):
 #             # find oxygens only bonded to current phosphorus
 #             bonded_oxygens_1 = [o for o in bonded_oxygens if len(o.get_bonded_heavy_atoms())==1]
 #             # find the closest oxygen ...
-#             closest_oxygen = min(bonded_oxygens_1, 
+#             closest_oxygen = min(bonded_oxygens_1,
-#                                  key= lambda o:Source.calculations.squared_distance(atom,o))
+#                                  key= lambda o:propka.calculations.squared_distance(atom,o))
 #             # ... and set it to O.2
 #             set_type(closest_oxygen,'O.2')
 
@@ -378,20 +383,20 @@ def are_atoms_planar(atoms):
     v1 = vector(atom1=atoms[0], atom2=atoms[1])
     v2 = vector(atom1=atoms[0], atom2=atoms[2])
     n = (v1**v2).rescale(1.0)
-    
+
     margin = 0.20
     for b in atoms[3:]:
         v = vector(atom1=atoms[0], atom2=b).rescale(1.0)
         #print(atoms[0],abs(v*n) )
         if abs(v*n)>margin:
             return False
-            
+
     return True
 
 def is_aromatic_ring(atoms):
     if len(atoms)<5:
         return False
-    
+
     for i in range(len(atoms)):
         if not are_atoms_planar(atoms[i:]+atoms[:i]):
             return False

propka/ligand_pka_values.py

@@ -1,6 +1,9 @@
 #!/usr/bin/env python
 
-import Source.molecular_container, Source.calculations, Source.calculations, Source.parameters, Source.pdb, Source.lib, os, subprocess, sys
+from __future__ import division
+from __future__ import print_function
+
+import propka.molecular_container, propka.calculations, propka.calculations, propka.parameters, propka.pdb, propka.lib, os, subprocess, sys
 
 class ligand_pka_values:
     def __init__(self, parameters):
@@ -14,14 +17,14 @@ class ligand_pka_values:
         print(self.molconvert)
 
         return
-    
+
-    
+
     def find_in_path(self, program):
         path = os.environ.get('PATH').split(os.pathsep)
 
         l = [i for i in filter(lambda loc: os.access(loc, os.F_OK),
                                map(lambda dir: os.path.join(dir, program),path))]
-    
+
         if len(l) == 0:
             print('Error: Could not find %s. Please make sure that it is found in the path.'%program)
             sys.exit(-1)
@@ -29,14 +32,14 @@ class ligand_pka_values:
         return l[0]
 
     def get_marvin_pkas_for_pdb_file(self, file, no_pkas=10, min_pH =-10, max_pH=20):
-        molecule = Source.molecular_container.Molecular_container(file)
+        molecule = propka.molecular_container.Molecular_container(file)
         self.get_marvin_pkas_for_molecular_container(molecule, no_pkas=no_pkas, min_pH =min_pH, max_pH=max_pH)
         return
 
     def get_marvin_pkas_for_molecular_container(self, molecule, no_pkas=10, min_pH =-10, max_pH=20):
         for name in molecule.conformation_names:
             filename = '%s_%s'%(molecule.name,name)
-            self.get_marvin_pkas_for_conformation_container(molecule.conformations[name], name=filename, reuse=molecule.options.reuse_ligand_mol2_file, 
+            self.get_marvin_pkas_for_conformation_container(molecule.conformations[name], name=filename, reuse=molecule.options.reuse_ligand_mol2_file,
                                                             no_pkas=no_pkas, min_pH =min_pH, max_pH=max_pH)
 
         return
@@ -50,7 +53,7 @@ class ligand_pka_values:
 
     def get_marvin_pkas_for_atoms(self, atoms, name='temp', reuse=False, no_pkas=10, min_pH =-10, max_pH=20):
         # do one molecule at the time so we don't confuse marvin
-        molecules = Source.lib.split_atoms_into_molecules(atoms)
+        molecules = propka.lib.split_atoms_into_molecules(atoms)
         for i in range(len(molecules)):
             filename = '%s_%d.mol2'%(name, i+1)
             self.get_marvin_pkas_for_molecule(molecules[i], filename=filename, reuse=reuse, no_pkas=no_pkas, min_pH =min_pH, max_pH=max_pH)
@@ -61,18 +64,18 @@ class ligand_pka_values:
     def get_marvin_pkas_for_molecule(self, atoms, filename='__tmp_ligand.mol2', reuse=False, no_pkas=10, min_pH =-10, max_pH=20):
         # print out structure unless we are using user-modified structure
         if not reuse:
-            Source.pdb.write_mol2_for_atoms(atoms, filename)
+            propka.pdb.write_mol2_for_atoms(atoms, filename)
         # check that we actually have a file to work with
         if not os.path.isfile(filename):
             print('Warning: Didn\'t find a user-modified file \'%s\' - generating one'%filename)
-            Source.pdb.write_mol2_for_atoms(atoms, filename)
+            propka.pdb.write_mol2_for_atoms(atoms, filename)
 
 
 
         # Marvin
         # calculate pKa values
         options = 'pka -a %d -b %d --min %f --max %f -d large'%(no_pkas, no_pkas, min_pH, max_pH)
-        (output,errors) = subprocess.Popen([self.cxcalc, filename]+options.split(), 
+        (output,errors) = subprocess.Popen([self.cxcalc, filename]+options.split(),
                                   stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate()
 
         if len(errors)>0:
@@ -86,7 +89,7 @@ class ligand_pka_values:
 
         # extract calculated pkas
         indices,pkas,types = self.extract_pkas(output)
-        
+
         # store calculated pka values
         for i in range(len(indices)):
             atoms[indices[i]].marvin_pka = pkas[i]
@@ -104,10 +107,10 @@ class ligand_pka_values:
         values = values.split('\t')
 
         # format values
-        types = [tags[i][0] for i in range(1,len(tags)-1) if len(values)>i and values[i] != ''] 
+        types = [tags[i][0] for i in range(1,len(tags)-1) if len(values)>i and values[i] != '']
         indices = [int(a)-1 for a in values[-1].split(',') if a !='']
         values = [float(v.replace(',','.')) for v in values[1:-1] if v != '']
-        
+
         if len(indices) != len(values) != len(types):
             raise Exception('Lengths of atoms and pka values mismatch')
 

propka/molecular_container.py

@@ -1,15 +1,19 @@
 #!/usr/bin/python
 #
-# Molecular container for storing all contents of pdb files 
+# Molecular container for storing all contents of pdb files
 #
 #
+from __future__ import division
+from __future__ import print_function
 
-import os, Source.pdb, sys, Source.version, Source.output, Source.conformation_container, Source.group, Source.lib
+import os, sys
+
+import propka.pdb, propka.version, propka.output, propka.conformation_container, propka.group, propka.lib
 
 class Molecular_container:
     def __init__(self, input_file, options=None):
         # printing out header before parsing input
-        Source.output.printHeader()
+        propka.output.printHeader()
 
         # set up some values
         self.options = options
@@ -21,11 +25,11 @@ class Molecular_container:
 
         # set the version
         if options:
-            parameters = Source.parameters.Parameters(self.options.parameters)
+            parameters = propka.parameters.Parameters(self.options.parameters)
         else:
-            parameters = Source.parameters.Parameters('propka.cfg')
+            parameters = propka.parameters.Parameters('propka.cfg')
         try:
-            exec('self.version = Source.version.%s(parameters)'%parameters.version)
+            exec('self.version = propka.version.%s(parameters)'%parameters.version)
         except:
             raise Exception('Error: Version %s does not exist'%parameters.version)
 
@@ -33,7 +37,7 @@ class Molecular_container:
         if input_file_extension[0:4] == '.pdb':
             # input is a pdb file
             # read in atoms and top up containers to make sure that all atoms are present in all conformations
-            [self.conformations, self.conformation_names] = Source.pdb.read_pdb(input_file, self.version.parameters,self)
+            [self.conformations, self.conformation_names] = propka.pdb.read_pdb(input_file, self.version.parameters,self)
             if len(self.conformations)==0:
                 print('Error: The pdb file does not seems to contain any molecular conformations')
                 sys.exit(-1)
@@ -41,7 +45,7 @@ class Molecular_container:
             self.top_up_conformations()
 
             # make a structure precheck
-            Source.pdb.protein_precheck(self.conformations, self.conformation_names)
+            propka.pdb.protein_precheck(self.conformations, self.conformation_names)
 
             # set up atom bonding and protonation
             self.version.setup_bonding_and_protonation(self)
@@ -58,12 +62,12 @@ class Molecular_container:
 
             # write out the input file
             filename = self.file.replace(input_file_extension,'.propka_input')
-            Source.pdb.write_input(self, filename)
+            propka.pdb.write_input(self, filename)
 
 
         elif input_file_extension == '.propka_input':
             #input is a propka_input file
-            [self.conformations, self.conformation_names] = Source.pdb.read_input(input_file, self.version.parameters, self)
+            [self.conformations, self.conformation_names] = propka.pdb.read_input(input_file, self.version.parameters, self)
 
             # Extract groups - this merely sets up the groups found in the input file
             self.extract_groups()
@@ -103,15 +107,15 @@ class Molecular_container:
         """ Identify the groups needed for pKa calculation """
         for name in self.conformation_names:
             self.conformations[name].extract_groups()
-            
+
         return
 
     def additional_setup_when_reading_input_file(self):
         for name in self.conformation_names:
             self.conformations[name].additional_setup_when_reading_input_file()
-            
+
         return
-        
+
 
     def calculate_pka(self):
         # calculate for each conformation
@@ -120,21 +124,21 @@ class Molecular_container:
 
         # find non-covalently coupled groups
         self.find_non_covalently_coupled_groups()
-        
+
         # find the average of the conformations
         self.average_of_conformations()
 
         # print out the conformation-average results
-        Source.output.printResult(self, 'AVR', self.version.parameters)
+        propka.output.printResult(self, 'AVR', self.version.parameters)
 
         return
 
     def average_of_conformations(self):
         # make a new configuration to hold the average values
-        avr_conformation = Source.conformation_container.Conformation_container(name='average',
+        avr_conformation = propka.conformation_container.Conformation_container(name='average',
-                                                                                parameters=self.conformations[self.conformation_names[0]].parameters, 
+                                                                                parameters=self.conformations[self.conformation_names[0]].parameters,
                                                                                 molecular_container=self)
-        
+
         for group in self.conformations[self.conformation_names[0]].get_titratable_groups_and_cysteine_bridges():
             # new group to hold average values
             avr_group = group.clone()
@@ -161,13 +165,13 @@ class Molecular_container:
 
     def write_pka(self, filename=None, reference="neutral", direction="folding", options=None):
         #for name in self.conformation_names:
-        #    Source.output.writePKA(self, self.version.parameters, filename='%s_3.1_%s.pka'%(self.name, name), 
+        #    propka.output.writePKA(self, self.version.parameters, filename='%s_3.1_%s.pka'%(self.name, name),
-        #                           conformation=name,reference=reference, 
+        #                           conformation=name,reference=reference,
         #                           direction=direction, options=options)
 
         # write out the average conformation
         filename=os.path.join('%s.pka'%(self.name))
-        
+
         # if the display_coupled_residues option is true,
         # write the results out to an alternative pka file
         if self.options.display_coupled_residues:
@@ -176,8 +180,8 @@ class Molecular_container:
         if hasattr(self.version.parameters, 'output_file_tag') and len(self.version.parameters.output_file_tag)>0:
             filename=os.path.join('%s_%s.pka'%(self.name,self.version.parameters.output_file_tag))
 
-        Source.output.writePKA(self, self.version.parameters, filename=filename,
+        propka.output.writePKA(self, self.version.parameters, filename=filename,
-                               conformation='AVR',reference=reference, 
+                               conformation='AVR',reference=reference,
                                direction=direction, options=options)
 
         return
@@ -185,7 +189,7 @@ class Molecular_container:
     def getFoldingProfile(self, conformation='AVR',reference="neutral", direction="folding", grid=[0., 14., 0.1], options=None):
         # calculate stability profile
         profile = []
-        for ph in Source.lib.make_grid(*grid):
+        for ph in propka.lib.make_grid(*grid):
             ddg = self.conformations[conformation].calculate_folding_energy( pH=ph, reference=reference)
             #print(ph,ddg)
             profile.append([ph, ddg])
@@ -196,29 +200,29 @@ class Molecular_container:
             opt = min(opt, point, key=lambda v:v[1])
 
         # find values within 80 % of optimum
-        range_80pct = [None, None] 
+        range_80pct = [None, None]
         values_within_80pct = [p[0] for p in profile if p[1]< 0.8*opt[1]]
         if len(values_within_80pct)>0:
             range_80pct = [min(values_within_80pct), max(values_within_80pct)]
 
         # find stability range
         stability_range = [None, None]
-        stable_values = [p[0] for p in profile if p[1]< 0.0]        
+        stable_values = [p[0] for p in profile if p[1]< 0.0]
-        
+
         if len(stable_values)>0:
             stability_range = [min(stable_values), max(stable_values)]
-        
+
         return profile, opt, range_80pct, stability_range
 
 
     def getChargeProfile(self, conformation='AVR', grid=[0., 14., .1]):
         charge_profile = []
-        for ph in Source.lib.make_grid(*grid):
+        for ph in propka.lib.make_grid(*grid):
             q_unfolded, q_folded = self.conformations[conformation].calculate_charge(self.version.parameters, pH=ph)
             charge_profile.append([ph, q_unfolded, q_folded])
 
         return charge_profile
-        
+
     def getPI(self, conformation='AVR', grid=[0., 14., 1], iteration=0):
         #print('staring',grid, iteration)
         # search
@@ -240,7 +244,7 @@ class Molecular_container:
 
         return pi_folded_value, pi_unfolded_value
 
-    
+
 
 if __name__ == '__main__':
     input_file = sys.argv[1]

propka/output.py

@@ -1,5 +1,10 @@
+
+from __future__ import division
+from __future__ import print_function
+
 import sys
-import Source.lib
+
+import propka.lib
 
 
 def printHeader():
@@ -17,7 +22,7 @@ def writePDB(protein, file=None, filename=None, include_hydrogens=False, options
     """
     Write the residue to the new pdbfile
     """
-    
+
     if file == None:
       # opening file if not given
       if filename == None:
@@ -124,7 +129,7 @@ def getDeterminantSection(protein, conformation, parameters):
     # printing determinants
     for chain in protein.conformations[conformation].chains:
         for residue_type in parameters.write_out_order:
-            groups = [g for g in protein.conformations[conformation].groups if g.atom.chainID == chain] 
+            groups = [g for g in protein.conformations[conformation].groups if g.atom.chainID == chain]
             for group in groups:
                 if group.residue_type == residue_type:
                     str += "%s" % ( group.getDeterminantString(parameters.remove_penalised_group) )
@@ -158,8 +163,8 @@ def getFoldingProfileSection(protein, conformation='AVR', direction="folding", r
     str += "\n"
     str += "Free energy of %9s (kcal/mol) as a function of pH (using %s reference)\n" % (direction, reference)
 
-    profile, [pH_opt, dG_opt], [dG_min, dG_max], [pH_min, pH_max] = protein.getFoldingProfile(conformation=conformation, 
+    profile, [pH_opt, dG_opt], [dG_min, dG_max], [pH_min, pH_max] = protein.getFoldingProfile(conformation=conformation,
-                                                                                 reference=reference, 
+                                                                                 reference=reference,
                                                                                  direction=direction, grid=[0., 14., 0.1], options=options)
     if profile == None:
       str += "Could not determine folding profile\n"
@@ -183,7 +188,7 @@ def getFoldingProfileSection(protein, conformation='AVR', direction="folding", r
       str += "Could not determine the pH-range where the free energy is negative\n\n"
     else:
       str += "The free energy is negative in the range %4.1lf - %4.1lf\n\n" % (pH_min, pH_max)
-     
+
 
     return str
 
@@ -209,7 +214,7 @@ def getChargeProfileSection(protein, conformation='AVR', options=None):
       str += "Could not determine the pI\n\n"
     else:
       str += "The pI is %5.2lf (folded) and %5.2lf (unfolded)\n" % (pI_pro, pI_mod)
-     
+
 
     return str
 
@@ -300,8 +305,8 @@ def getReferencesHeader():
     str += "   Journal of Chemical Theory and Computation, 7(2):525-537 (2011)\n"
     str += "   \n"
     str += "   Improved Treatment of Ligands and Coupling Effects in Empirical Calculation\n"
-    str += "    and Rationalization of pKa Values\n" 
+    str += "    and Rationalization of pKa Values\n"
-    str += "   Chresten R. Sondergaard, Mats H.M. Olsson, Michal Rostkowski, and Jan H. Jensen\n" 
+    str += "   Chresten R. Sondergaard, Mats H.M. Olsson, Michal Rostkowski, and Jan H. Jensen\n"
     str += "   Journal of Chemical Theory and Computation, (2011)\n"
     str += "   \n"
     str += "-------------------------------------------------------------------------------------------------------\n"
@@ -357,13 +362,13 @@ def getTheLine():
 
 # Interaction maps
 def make_interaction_map(name, list, interaction):
-    """ Print out an interaction map named 'name' of the groups in 'list' 
+    """ Print out an interaction map named 'name' of the groups in 'list'
     based on the function 'interaction' """
-    
+
     # return an empty string, if the list is empty
     if len(list)==0:
         return ''
-    
+
     # for long list, use condensed formatting
     if len(list)>10:
         res = 'Condensed form:\n'
@@ -377,7 +382,7 @@ def make_interaction_map(name, list, interaction):
     res = '%s\n%12s'%(name,'')
     for g in list:
         res += '%9s | '%g.label
-    
+
     # do the map
     for g1 in list:
         res += '\n%-12s'%(g1.label)
@@ -386,6 +391,6 @@ def make_interaction_map(name, list, interaction):
             if interaction(g1, g2):
                 tag = '    X     '
             res += '%10s| '%tag
-    
+
     return res
-        
+

propka/parameters.py

@@ -1,7 +1,12 @@
+
+from __future__ import division
+from __future__ import print_function
+
 import math
-import Source.lib as lib
+import propka.lib as lib
 import sys, os
 
+import pkg_resources
 
 # names and types of all key words in configuration file
 matrices =            ['interaction_matrix']
@@ -45,16 +50,15 @@ class Parameters:
         #self.print_interaction_parameters_latex()
         #####self.print_interactions_latex()
         #sys.exit(0)
-        
+
-        
+
         return
 
 
     def read_parameters(self, file):
         # try to locate the parameters file
         try:
-            path = os.path.dirname(__file__)
+            ifile = pkg_resources.resource_filename(__name__, file)
-            ifile = os.path.join(path,'../'+file)
             input = lib.open_file_for_reading(ifile)
         except:
             input = lib.open_file_for_reading(file)
@@ -94,7 +98,7 @@ class Parameters:
         elif len(words)==3 and words[0] in string_dictionaries:
             self.parse_to_string_dictionary(words)
 
-            
+
         #print(words)
 
         return
@@ -226,25 +230,25 @@ class Parameters:
         print("""
 Titratable:
 CG  ARG
-C2N ARG 
+C2N ARG
 N30 N+/LYS
 N31 N+/LYS
 N32 N+/LYS
-N33 N+/LYS 
+N33 N+/LYS
 NAR HIS
 OCO COO
 OP  TYR/SER?
-SH  CYS 
+SH  CYS
 
 Non-titratable:
 NP1 AMD?
 OH  ROH
 O3  ?
-CL  
+CL
-F   
+F
-NAM 
+NAM
-N1  
+N1
-O2 
+O2
 """)
         return
 
@@ -283,7 +287,7 @@ O2
             'N1' :[],
             'O2' :[]}
 
-   
+
         s = """
 \\begin{longtable}{lllll}
 \\caption{Ligand interaction parameters. For interactions not listed, the default value of %s is applied.}
@@ -296,11 +300,11 @@ Group1 & Group2 & Interaction & c1 &c2 \\\\
 
 \\multicolumn{5}{l}{\\emph{continued from the previous page}}\\\\
 \\toprule
-Group1 & Group2 & Interaction & c1 &c2 \\\\ 
+Group1 & Group2 & Interaction & c1 &c2 \\\\
 \\midrule
 \\endhead
 
-\\midrule 
+\\midrule
 \\multicolumn{5}{r}{\\emph{continued on the next page}}\\\\
 \\endfoot
 
@@ -332,7 +336,7 @@ Group1 & Group2 & Interaction & c1 &c2 \\\\
         agroups = ['COO', 'HIS', 'CYS', 'TYR', 'SER', 'N+', 'LYS', 'AMD', 'ARG', 'TRP', 'ROH', 'CG', 'C2N', 'N30', 'N31', 'N32', 'N33', 'NAR', 'OCO', 'NP1', 'OH', 'O3', 'CL', 'F', 'NAM', 'N1', 'O2', 'OP', 'SH']
         lgroups = ['CG', 'C2N', 'N30', 'N31', 'N32', 'N33', 'NAR', 'OCO', 'NP1', 'OH', 'O3', 'CL', 'F', 'NAM', 'N1', 'O2', 'OP', 'SH']
 
-   
+
         s = """
 \\begin{longtable}{%s}
 \\caption{Ligand interaction parameters. For interactions not listed, the default value of %s is applied.}
@@ -345,11 +349,11 @@ Group1 & Group2 & Interaction & c1 &c2 \\\\
 
 \\multicolumn{5}{l}{\\emph{continued from the previous page}}\\\\
 \\toprule
-Group1 & Group2 & Interaction & c1 &c2 \\\\ 
+Group1 & Group2 & Interaction & c1 &c2 \\\\
 \\midrule
 \\endhead
 
-\\midrule 
+\\midrule
 \\multicolumn{5}{r}{\\emph{continued on the next page}}\\\\
 \\endfoot
 
@@ -399,14 +403,14 @@ class Interaction_matrix:
                 self.dictionary[group][new_group] = self.value
                 self.dictionary[new_group][group] = self.value
 
-        
+
         return
 
     def get_value(self, item1, item2):
         try:
             return self.dictionary[item1][item2]
         except:
-            return None 
+            return None
 
     def __getitem__(self, group):
         if group not in self.dictionary.keys():
@@ -430,7 +434,7 @@ class Interaction_matrix:
 
         return s
 #         ks = ['COO', 'SER', 'ARG', 'LYS', 'HIS', 'AMD', 'CYS', 'TRP','ROH','TYR','N+','CG', 'C2N', 'N30', 'N31', 'N32', 'N33', 'NAR', 'OCO', 'NP1', 'OH', 'O3', 'CL', 'F', 'NAM', 'N1', 'O2', 'OP', 'SH']
- 
+
 #         p = ''
 #         n=0
 #         for i in range(len(ks)):
@@ -442,7 +446,7 @@ class Interaction_matrix:
 
 #         print('total',n,len(ks))
 #         return p
-    
+
 
 
 class Pair_wise_matrix:
@@ -451,7 +455,7 @@ class Pair_wise_matrix:
         self.dictionary = {}
         self.default = [0.0, 0.0]
         return
-    
+
     def add(self,words):
         # assign the default value
         if len(words)==3 and words[0]=='default':
@@ -462,7 +466,7 @@ class Pair_wise_matrix:
         g1 = words[0]
         g2 = words[1]
         v = [float(words[2]), float(words[3])]
-        
+
         self.insert(g1,g2,v)
         self.insert(g2,g1,v)
 
@@ -482,7 +486,7 @@ class Pair_wise_matrix:
         return
 
     def get_value(self, item1, item2):
-        
+
         try:
             return self.dictionary[item1][item2]
         except:
@@ -504,7 +508,7 @@ class Pair_wise_matrix:
                 s += '%s %s %s\n'%(k1,k2,self[k1][k2])
 
         return s
-    
+
 
 
 

propka/pdb.py

@@ -1,6 +1,12 @@
-import string, sys, copy, Source.lib
+
-from Source.atom import Atom
+from __future__ import division
-from Source.conformation_container import Conformation_container
+from __future__ import print_function
+
+import string, sys, copy
+
+import propka.lib
+from propka.atom import Atom
+from propka.conformation_container import Conformation_container
 
 expected_atom_numbers = {'ALA':5,
                          'ARG':11,
@@ -35,23 +41,23 @@ def read_pdb(pdb_file, parameters, molecule):
         conformations[name].add_atom(atom)
 
     # make a sorted list of conformation names
-    names = sorted(conformations.keys(), key=Source.lib.conformation_sorter)
+    names = sorted(conformations.keys(), key=propka.lib.conformation_sorter)
- 
+
     return [conformations, names]
 
 def protein_precheck(conformations, names):
-    
+
     for name in names:
         atoms = conformations[name].atoms
-        
+
         res_ids = []
         [res_ids.append(resid_from_atom(a)) for a in atoms if not res_ids.count(resid_from_atom(a))]
-        
+
         for res_id in res_ids:
             res_atoms = [a for a in atoms if resid_from_atom(a) == res_id and a.element != 'H']
             resname = res_atoms[0].resName
             residue_label = '%3s%5s'%(resname, res_id)
-            
+
             # ignore ligand residues
             if resname not in expected_atom_numbers:
                 continue
@@ -65,7 +71,7 @@ def protein_precheck(conformations, names):
             # check number of atoms in residue
             if len(res_atoms) != expected_atom_numbers[resname]:
                 print('Warning: Unexpected number (%d) of atoms in residue %s in conformation %s'%(len(res_atoms),residue_label, name))
-            
+
     return
 
 def resid_from_atom(a):
@@ -74,7 +80,7 @@ def resid_from_atom(a):
 
 def get_atom_lines_from_pdb(pdb_file, ignore_residues = [], keep_protons=False, tags = ['ATOM  ', 'HETATM'], chains=None):
 
-    lines = Source.lib.open_file_for_reading(pdb_file).readlines()
+    lines = propka.lib.open_file_for_reading(pdb_file).readlines()
     nterm_residue = 'next_residue'
     old_residue = None
     terminal = None
@@ -88,7 +94,7 @@ def get_atom_lines_from_pdb(pdb_file, ignore_residues = [], keep_protons=False,
         if tag == 'MODEL ':
             model = int(line[6:])
             nterm_residue = 'next_residue'
-        
+
         if tag == 'TER   ':
             nterm_residue = 'next_residue'
 
@@ -102,21 +108,21 @@ def get_atom_lines_from_pdb(pdb_file, ignore_residues = [], keep_protons=False,
                 continue
             if chains and line[21] not in chains:
                 continue
-            
+
-            # set the Nterm residue number - nessecary because we may need to 
+            # set the Nterm residue number - nessecary because we may need to
             # identify more than one N+ group for structures with alt_conf tags
             if nterm_residue == 'next_residue' and tag == 'ATOM  ':
-                # make sure that we reached a new residue - nessecary if OXT is not the last atom in 
+                # make sure that we reached a new residue - nessecary if OXT is not the last atom in
-                # the previous residue 
+                # the previous residue
                 if old_residue != residue_number:
                     nterm_residue = residue_number
                     old_residue = None
 
 
-            # Identify the configuration 
+            # Identify the configuration
             # convert digits to letters
             if alt_conf_tag in '123456789':
-                alt_conf_tag = chr(ord(alt_conf_tag)+16) 
+                alt_conf_tag = chr(ord(alt_conf_tag)+16)
             if alt_conf_tag == ' ':
                 alt_conf_tag = 'A'
             conformation = '%d%s'%(model, alt_conf_tag)
@@ -145,10 +151,10 @@ def get_atom_lines_from_pdb(pdb_file, ignore_residues = [], keep_protons=False,
 def write_pdb(conformation, filename):
     write_pdb_for_atoms(conformation.atoms, filename)
     return
-    
+
 def write_pdb_for_atoms(atoms, filename, make_conect_section=False):
-    out = Source.lib.open_file_for_writing(filename)
+    out = propka.lib.open_file_for_writing(filename)
-    
+
     for atom in atoms:
         out.write(atom.make_pdb_line())
 
@@ -171,7 +177,7 @@ def write_mol2_for_atoms(atoms, filename):
     for i in range(len(atoms)):
         atoms_section += atoms[i].make_mol2_line(i+1)
 
-        
+
     bonds_section = '@<TRIPOS>BOND\n'
     id = 1
     for i in range(len(atoms)):
@@ -180,12 +186,12 @@ def write_mol2_for_atoms(atoms, filename):
                 type = get_bond_order(atoms[i],atoms[j])
                 bonds_section += '%7d %7d %7d %7s\n'%(id, i+1, j+1, type)
                 id+=1
-    
+
     substructure_section = '@<TRIPOS>SUBSTRUCTURE\n\n'
     if len(atoms)>0:
         substructure_section = '@<TRIPOS>SUBSTRUCTURE\n%-7d %10s %7d\n'%(atoms[0].resNumb,atoms[0].resName,atoms[0].numb)
 
-    out = Source.lib.open_file_for_writing(filename)
+    out = propka.lib.open_file_for_writing(filename)
     out.write(header%(len(atoms),id-1))
     out.write(atoms_section)
     out.write(bonds_section)
@@ -209,14 +215,14 @@ def get_bond_order(atom1, atom2):
 
     if '.ar' in atom1.sybyl_type and '.ar' in atom2.sybyl_type:
         type = 'ar'
-          
+
 
     return type
 
 
 
 def write_input(molecular_container, filename):
-    out = Source.lib.open_file_for_writing(filename)
+    out = propka.lib.open_file_for_writing(filename)
 
     for conformation_name in molecular_container.conformation_names:
         out.write('MODEL %s\n'%conformation_name)
@@ -250,14 +256,14 @@ def read_input(input_file, parameters,molecule):
         conformations[name].add_atom(atom)
 
     # make a sorted list of conformation names
-    names = sorted(conformations.keys(), key=Source.lib.conformation_sorter)
+    names = sorted(conformations.keys(), key=propka.lib.conformation_sorter)
- 
+
     return [conformations, names]
 
 
 
 def get_atom_lines_from_input(input_file, tags = ['ATOM  ','HETATM']):
-    lines = Source.lib.open_file_for_reading(input_file).readlines()
+    lines = propka.lib.open_file_for_reading(input_file).readlines()
     conformation = ''
 
     atoms = {}
@@ -278,7 +284,7 @@ def get_atom_lines_from_input(input_file, tags = ['ATOM  ','HETATM']):
             atom.is_protonated = True
             atoms[atom.numb] = atom
             numbers.append(atom.numb)
-            
+
         # found bonding information - apply it
         if tag == 'CONECT' and len(line)>14:
             conect_numbers = [line[i:i+5] for i in range(6, len(line)-1, 5)]
@@ -301,7 +307,7 @@ def get_atom_lines_from_input(input_file, tags = ['ATOM  ','HETATM']):
             center_atom = atoms[int(conect_numbers[0])]
             for n in conect_numbers[1:]:
                 cg = atoms[int(n)]
-                center_atom.group.couple_covalently(cg.group)            
+                center_atom.group.couple_covalently(cg.group)
 
         # found info on non-covalent coupling
         if tag == 'NCOUPL' and len(line)>14:
@@ -309,8 +315,8 @@ def get_atom_lines_from_input(input_file, tags = ['ATOM  ','HETATM']):
             center_atom = atoms[int(conect_numbers[0])]
             for n in conect_numbers[1:]:
                 cg = atoms[int(n)]
-                center_atom.group.couple_non_covalently(cg.group)            
+                center_atom.group.couple_non_covalently(cg.group)
-            
+
         # this conformation is done - yield the atoms
         if tag == 'ENDMDL':
             for n in numbers:

propka/propka.cfg

@@ -0,0 +1,398 @@
+# PropKa configuration file
+
+version version_A
+
+# Model pKa values
+model_pkas C-   3.20
+model_pkas ASP  3.80
+model_pkas GLU  4.50
+model_pkas HIS  6.50
+model_pkas CYS  9.00
+model_pkas TYR 10.00
+model_pkas LYS 10.50
+model_pkas ARG 12.50
+#model_pkas SER 14.20 Jack Kyte: Structure in Protein Chemistry, 1995, Garland Publishing, Inc New York and London
+model_pkas N+   8.00
+model_pkas CG  11.50
+model_pkas C2N 11.50
+model_pkas N30 10.00
+model_pkas N31 10.00
+model_pkas N32 10.00
+model_pkas N33 10.00
+model_pkas NAR  5.00
+model_pkas OCO  4.50
+model_pkas SH  10.00
+model_pkas OP   6.00
+
+# Custom ligand pKa values
+#  P. Acharya, P. Cheruku, S. Chatterjee, S. Acharya, and, J. Chattopadhyaya:
+#  Measurement of Nucleobase pKa Values in Model Mononucleotides 
+#  Shows RNA-RNA Duplexes To Be More Stable than DNA-DNA Duplexes
+#  Journal of the American Chemical Society 2004 126 (9), 2862-2869
+# 
+custom_model_pkas DA-N1   3.82
+custom_model_pkas DA-N3   3.82
+custom_model_pkas DA-N7   3.82
+custom_model_pkas DA-OP1  1.00
+custom_model_pkas DA-OP2  1.00
+
+custom_model_pkas DG-N1   9.59
+custom_model_pkas DG-N3   9.59
+custom_model_pkas DG-N7   9.59
+custom_model_pkas DG-OP1  1.00
+custom_model_pkas DG-OP2  1.00
+
+custom_model_pkas DC-N3   4.34
+custom_model_pkas DC-OP1  1.00
+custom_model_pkas DC-OP2  1.00
+
+custom_model_pkas DT-N3  10.12
+custom_model_pkas DT-OP1  1.00
+custom_model_pkas DT-OP2  1.00
+
+
+# protein group mapping
+protein_group_mapping ASP-CG  COO
+protein_group_mapping GLU-CD  COO
+protein_group_mapping HIS-CG  HIS
+protein_group_mapping CYS-SG  CYS
+protein_group_mapping TYR-OH  TYR
+protein_group_mapping LYS-NZ  LYS
+protein_group_mapping ARG-CZ  ARG
+#protein_group_mapping SER-OG  SER
+protein_group_mapping THR-OG1 ROH
+protein_group_mapping SER-OG  ROH#
+protein_group_mapping ASN-CG  AMD
+protein_group_mapping GLN-CD  AMD
+protein_group_mapping TRP-NE1 TRP
+
+
+# matrix for propka interactions
+#   'N'   non-iterative interaction
+#   'I'   iterative interaction
+#   '-'   no interaction
+                      #CYS
+interaction_matrix CYS I#N+ 
+interaction_matrix N+  N I#HIS 
+interaction_matrix HIS I N I#LYS 
+interaction_matrix LYS N N N I#AMD 
+interaction_matrix AMD N - N - -#COO 
+interaction_matrix COO I N I N N I#ARG 
+interaction_matrix ARG N N N N - N I#TRP 
+interaction_matrix TRP N - - - - N - -#ROH 
+interaction_matrix ROH N - - - - N - - -#TYR 
+interaction_matrix TYR N I I I N N N N N I#SER
+interaction_matrix SER N N N N N N I N N N I #CG
+interaction_matrix CG  N N N N - N I - - N I I#C2N
+interaction_matrix C2N N N N N - N I - - N I I I#N30
+interaction_matrix N30 N I N N - N N - - I N I I I#N31
+interaction_matrix N31 N I N N - N N - - I N I I I I#N32
+interaction_matrix N32 N I N N - N N - - I N I I I I I#N33
+interaction_matrix N33 N I N N - N N - - I N I I I I I I#NAR
+interaction_matrix NAR I N I I N I N - - I N N N N N N N I#OCO
+interaction_matrix OCO I N I N N I N N N N N N N N N N N I I#NP1
+interaction_matrix NP1 N - N - - N - - - N N - - - - - - N N -#OH
+interaction_matrix OH  N - - - - N - - - N N - - - - - - - N - -#O3
+interaction_matrix O3  N - N - - N - - - N N - - - - - - N N - - -#CL
+interaction_matrix CL  N - N - - N - - - N N - - - - - - N N - - - -#F
+interaction_matrix F   N - N - - N - - - N N - - - - - - N N - - - - -#NAM
+interaction_matrix NAM N - N - - N - - - N N - - - - - - N N - - - - - -#N1
+interaction_matrix N1  N - N - - N - - - N N - - - - - - N N - - - - - - -#O2
+interaction_matrix O2  N - N - - N - - - N N - - - - - - N N - - - - - - - -#OP
+interaction_matrix OP  I N I N N I N N N N N N N N N N N I I N N N N N N N N I#SH
+interaction_matrix SH  I N N N N N N N N N N I I I I I I N N N N N N N N N N N I
+
+# Cutoff values for side chain interactions
+# default value
+sidechain_cutoffs default 3.0  4.0  
+# COO
+sidechain_cutoffs COO COO 2.5  3.5
+Sidechain_cutoffs COO SER 2.65 3.65
+sidechain_cutoffs COO ARG 1.85 2.85
+sidechain_cutoffs COO LYS 2.85 3.85
+sidechain_cutoffs COO HIS 2.0  3.0
+sidechain_cutoffs COO AMD 2.0  3.0
+sidechain_cutoffs COO TRP 2.0  3.0
+sidechain_cutoffs COO ROH 2.65 3.65
+sidechain_cutoffs COO TYR 2.65 3.65
+sidechain_cutoffs COO  N+ 2.85 3.85
+sidechain_cutoffs COO  CG 1.85 2.85
+sidechain_cutoffs COO C2N 1.85 2.85
+sidechain_cutoffs COO N30 2.85 3.85
+sidechain_cutoffs COO N31 2.85 3.85
+sidechain_cutoffs COO N32 2.85 3.85
+sidechain_cutoffs COO N33 2.85 3.85
+sidechain_cutoffs COO NAR 2.0  3.0
+sidechain_cutoffs COO OCO 2.5  3.5
+sidechain_cutoffs COO  OH 2.65 3.65
+sidechain_cutoffs COO NAM 2.0  3.0
+# SER
+sidechain_cutoffs SER SER 3.5  4.5
+sidechain_cutoffs SER ARG 2.5  4.0
+sidechain_cutoffs SER HIS 2.0  3.0
+sidechain_cutoffs SER AMD 2.5  3.5
+sidechain_cutoffs SER CYS 3.5  4.5
+sidechain_cutoffs SER TRP 2.5  3.5
+sidechain_cutoffs SER ROH 3.5  4.5
+sidechain_cutoffs SER  CG 2.5  4.0
+sidechain_cutoffs SER C2N 2.5  4.0
+sidechain_cutoffs SER NAR 2.0  3.0
+sidechain_cutoffs SER  OH 3.5  4.5
+sidechain_cutoffs SER  SH 3.5  4.5
+sidechain_cutoffs SER TYR 3.5  4.5
+sidechain_cutoffs SER  N+ 3.0  4.5
+sidechain_cutoffs SER NAM 2.5  3.5
+# ARG
+sidechain_cutoffs ARG CYS 2.5  4.0
+sidechain_cutoffs ARG TYR 2.5  4.0
+sidechain_cutoffs ARG OCO 1.85 2.85
+sidechain_cutoffs ARG  SH 2.5  4.0
+# HIS
+sidechain_cutoffs HIS AMD 2.0  3.0
+sidechain_cutoffs HIS TYR 2.0  3.0
+sidechain_cutoffs HIS OCO 2.0  3.0
+# CYS
+sidechain_cutoffs CYS CYS 3.0  5.0
+sidechain_cutoffs CYS TRP 2.5  3.5
+sidechain_cutoffs CYS ROH 3.5  4.5
+sidechain_cutoffs CYS AMD 2.5  3.5
+sidechain_cutoffs CYS TYR 3.5  4.5
+sidechain_cutoffs CYS  N+ 3.0  4.5
+sidechain_cutoffs CYS  CG 2.5  4.0
+sidechain_cutoffs CYS C2N 2.5  4.0
+sidechain_cutoffs CYS N30 3.0  4.5
+sidechain_cutoffs CYS N31 3.0  4.5
+sidechain_cutoffs CYS N32 3.0  4.5
+sidechain_cutoffs CYS N33 3.0  4.5
+sidechain_cutoffs CYS  OH 3.5  4.5
+sidechain_cutoffs CYS NAM 2.5  3.5
+sidechain_cutoffs CYS  SH 3.0  5.0
+# TYR
+sidechain_cutoffs TYR TYR 3.5  4.5
+sidechain_cutoffs TYR  N+ 3.0  4.5
+sidechain_cutoffs TYR AMD 2.5  3.5
+sidechain_cutoffs TYR TRP 2.5  3.5
+sidechain_cutoffs TYR ROH 3.5  4.5
+sidechain_cutoffs TYR  CG 2.5  4.0
+sidechain_cutoffs TYR C2N 2.5  4.0
+sidechain_cutoffs TYR OCO 2.65 3.65
+sidechain_cutoffs TYR NAR 2.0  3.0
+sidechain_cutoffs TYR  OH 3.5  4.5
+sidechain_cutoffs TYR NAM 2.5  3.5
+sidechain_cutoffs TYR  SH 3.5  4.5
+# N+
+sidechain_cutoffs  N+ OCO 2.85 3.85
+sidechain_cutoffs  N+  SH 3.0  4.5
+# LYS
+sidechain_cutoffs LYS OCO 2.85 3.85
+# OCO
+sidechain_cutoffs OCO OCO 2.5  3.5
+sidechain_cutoffs OCO TRP 2.0  3.0
+sidechain_cutoffs OCO ROH 2.65 3.65
+sidechain_cutoffs OCO AMD 2.0  3.0
+sidechain_cutoffs OCO  CG 1.85 2.85
+sidechain_cutoffs OCO C2N 1.85 2.85
+sidechain_cutoffs OCO N30 2.85 3.85
+sidechain_cutoffs OCO N31 2.85 3.85
+sidechain_cutoffs OCO N32 2.85 3.85
+sidechain_cutoffs OCO N33 2.85 3.85
+sidechain_cutoffs OCO NAR 2.0  3.0 
+sidechain_cutoffs OCO  OH 2.65 3.65
+sidechain_cutoffs OCO NAM 2.0  3.0
+# NAR
+sidechain_cutoffs NAR AMD 2.0  3.0
+# SH
+sidechain_cutoffs  SH ROH 3.5  4.5
+sidechain_cutoffs  SH TRP 2.5  3.5
+sidechain_cutoffs  SH AMD 2.5  3.5
+sidechain_cutoffs  SH NAM 2.5  3.5
+sidechain_cutoffs  SH  CG 2.5  4.0
+sidechain_cutoffs  SH C2N 2.5  4.0
+sidechain_cutoffs  SH  OH 3.5  4.5
+sidechain_cutoffs  SH  SH 3.0  5.0
+
+
+
+# Maximal interaction energies for side chains
+sidechain_interaction 0.85  
+
+# Angular dependent sidechain interactions
+angular_dependent_sidechain_interactions HIS
+angular_dependent_sidechain_interactions ARG
+angular_dependent_sidechain_interactions AMD
+angular_dependent_sidechain_interactions TRP
+
+# exception interaction values
+COO_HIS_exception 1.60
+OCO_HIS_exception 1.60
+CYS_HIS_exception 1.60
+CYS_CYS_exception 3.60
+
+# Coulomb interaction parameters
+coulomb_cutoff1  4.0  
+coulomb_cutoff2 10.0
+coulomb_diel    80.0
+
+# Backbone hydrogen bond parameters 
+backbone_NH_hydrogen_bond COO -0.85 2.00 3.00
+#backbone_NH_hydrogen_bond C-  -0.85 2.00 3.00
+backbone_NH_hydrogen_bond CYS -0.85 3.00 4.00
+backbone_NH_hydrogen_bond TYR -0.85 2.20 3.20
+backbone_NH_hydrogen_bond OCO -0.85 2.00 3.50
+backbone_NH_hydrogen_bond NAR -0.85 2.00 3.50
+
+backbone_CO_hydrogen_bond HIS  0.85 2.00 3.00
+backbone_CO_hydrogen_bond OCO  0.85 3.00 4.00
+backbone_CO_hydrogen_bond CG   0.85 2.00 4.00
+backbone_CO_hydrogen_bond C2N  0.85 2.00 4.00
+backbone_CO_hydrogen_bond N30  0.85 2.00 4.00
+backbone_CO_hydrogen_bond N31  0.85 2.00 4.00
+backbone_CO_hydrogen_bond N32  0.85 2.00 4.00
+backbone_CO_hydrogen_bond N33  0.85 2.00 4.00
+backbone_CO_hydrogen_bond NAR  0.85 2.00 3.50
+
+# Group charges
+charge COO -1
+charge HIS +1
+charge CYS -1
+charge TYR -1
+charge LYS +1
+charge ARG +1
+charge N+  +1
+charge C-  -1
+charge OCO -1
+charge SER -1
+charge CG  +1
+charge C2N +1
+charge N30 +1
+charge N31 +1
+charge N32 +1
+charge N33 +1
+charge NAR +1
+charge SH  -1
+charge OP  -1
+
+# list of acids
+acid_list ASP
+acid_list GLU
+acid_list CYS
+acid_list TYR
+acid_list SER
+acid_list C-
+acid_list OCO
+acid_list OP
+acid_list SH
+
+# list of bases
+base_list ARG
+base_list LYS
+base_list HIS
+base_list N+
+base_list CG
+base_list C2N
+base_list N30
+base_list N31
+base_list N32
+base_list N33
+base_list NAR
+
+# list of groups used in backbone reorganisation calculations
+backbone_reorganisation_list ASP
+backbone_reorganisation_list GLU
+
+# Residues that should be ignored
+ignore_residues HOH
+ignore_residues H2O
+ignore_residues HOH
+ignore_residues SO4
+ignore_residues PO4
+ignore_residues PEG
+ignore_residues EPE
+#ignore_residues NAG
+ignore_residues TRS
+
+# Relative Van der Waals volume parameters for the radial volume model
+# Radii adopted from Bondi, A. (1964). "Van der Waals Volumes and Radii". J. Phys. Chem. 68 (3): 441-51
+VanDerWaalsVolume C  1.40  # radius: 1.70, volume: 20.58 all 'C' and 'CA' atoms
+VanDerWaalsVolume C4 2.64  # 38.79 hydrodphobic carbon atoms + unidentified atoms
+VanDerWaalsVolume N  1.06  # radius: 1.55, volume: 15.60 all nitrogen atoms
+VanDerWaalsVolume O  1.00  # radius: 1.52, volume: 14.71 all oxygen atoms
+VanDerWaalsVolume S  1.66  # radius: 1.80, volume: 24.43 all sulphur atoms
+VanDerWaalsVolume F  0.90  # raidus: 1.47, volume: 13.30 for fluorine 
+VanDerWaalsVolume Cl 1.53  # radius: 1.75, volume: 22.44 for chlorine
+VanDerWaalsVolume P  1.66  # radius: 1.80, volume: 24.42 for phosphorus
+
+# Other desolvation parameters
+desolvationSurfaceScalingFactor	   0.25
+desolvationPrefactor		  -13.0
+desolvationAllowance		    0.0
+desolv_cutoff 			   20.0
+buried_cutoff 			   15.0
+Nmin	      			    280
+Nmax	      			    560
+
+# Ligand groups
+ligand_typing groups 
+min_bond_distance_for_hydrogen_bonds 4
+
+# covalent coupling
+coupling_max_number_of_bonds     3
+shared_determinants              0
+common_charge_centre             0
+remove_penalised_group	         1
+
+# non-covalent coupling
+max_intrinsic_pKa_diff         2.0
+min_interaction_energy         0.5
+max_free_energy_diff           1.0
+min_swap_pka_shift             1.0
+min_pka                        0.0
+max_pka                       10.0
+pH                        variable
+reference                  neutral
+
+# ions
+ions 1P   1 # generic charged atoms
+ions 2P   2
+ions 1N  -1
+ions 2N  -2
+
+ions MG   2 #Magnesium Ion
+ions CA   2 #Calcium Ion
+ions ZN   2 #Zinc Ion
+ions NA   1 #Sodium Ion
+ions CL  -1 #Chloride Ion
+ions MN   2 #Manganese (ii) Ion
+ions K    1 #Potassium Ion
+ions CD   2 #Cadmium Ion
+ions FE   3 #Fe (iii) Ion
+ions SR   2 #Strontium Ion
+ions CU   2 #Copper (ii) Ion
+ions IOD -1 #Iodide Ion
+ions HG   2 #Mercury (ii) Ion
+ions BR  -1 #Bromide Ion
+ions CO   2 #Cobalt (ii) Ion
+ions NI   2 #Nickel (ii) Ion
+ions FE2  2 #Fe (ii) Ion
+
+# write out order of residues
+write_out_order ASP 
+write_out_order GLU
+write_out_order C-
+write_out_order HIS
+write_out_order CYS
+write_out_order TYR
+write_out_order LYS
+write_out_order ARG
+write_out_order SER
+write_out_order N+
+write_out_order CG
+write_out_order C2N
+write_out_order N30
+write_out_order N31
+write_out_order N32
+write_out_order N33
+write_out_order NAR
+write_out_order OCO
+write_out_order SH
+write_out_order OP

propka/protonate.py

@@ -1,11 +1,14 @@
 #!/usr/bin/python
 
-from Source.vector_algebra import *
+from __future__ import division
-import Source.bonds, Source.pdb, Source.atom
+from __future__ import print_function
+
+from propka.vector_algebra import *
+import propka.bonds, propka.pdb, propka.atom
 
 class Protonate:
     """ Protonates atoms using VSEPR theory """
-    
+
     def __init__(self, verbose=False):
         self.verbose=verbose
 
@@ -56,9 +59,9 @@ class Protonate:
                                 'HIS-ND1':1.0,
                                 'LYS-NZ':1.0,
                                 'N+':1.0,
-                                'C-':-1.0} 
+                                'C-':-1.0}
+
 
-        
         self.sybyl_charges = {'N.pl3':+1,
                               'N.3':+1,
                               'N.4':+1,
@@ -92,17 +95,17 @@ class Protonate:
         self.display('----- Protonation started -----')
         # Remove all currently present hydrogen atoms
         self.remove_all_hydrogen_atoms(molecules)
-             
+
         # protonate all atoms
         for name in molecules.conformation_names:
             non_H_atoms = molecules.conformations[name].get_non_hydrogen_atoms()
 
             for atom in non_H_atoms:
                 self.protonate_atom(atom)
-            
+
             # fix hydrogen names
             #self.set_proton_names(non_H_atoms)
-                    
+
         return
 
 
@@ -110,7 +113,7 @@ class Protonate:
         for name in molecular_container.conformation_names:
             molecular_container.conformations[name].atoms = molecular_container.conformations[name].get_non_hydrogen_atoms()
         return
-    
+
 
     def set_charge(self, atom):
         # atom is a protein atom
@@ -148,11 +151,11 @@ class Protonate:
         for heavy_atom in heavy_atoms:
             i = 1
             for bonded in heavy_atom.bonded_atoms:
-                
+
                 if bonded.element == 'H':
                     bonded.name+='%d'%i
                     i+=1
-                        
+
 
         return
 
@@ -160,7 +163,7 @@ class Protonate:
     def set_number_of_protons_to_add(self, atom):
         self.display('*'*10)
         self.display('Setting number of protons to add for',atom)
-        atom.number_of_protons_to_add  = 8 
+        atom.number_of_protons_to_add  = 8
         self.display('                  %4d'%8)
         atom.number_of_protons_to_add -= self.valence_electrons[atom.element]
         self.display('Valence eletrons: %4d'%-self.valence_electrons[atom.element])
@@ -177,7 +180,7 @@ class Protonate:
         return
 
     def set_steric_number_and_lone_pairs(self, atom):
-        
+
         # If we already did this, there is no reason to do it again
         if atom.steric_number_and_lone_pairs_set:
             return
@@ -195,10 +198,10 @@ class Protonate:
 
 
         atom.steric_number = 0
-        
+
         self.display('%65s: %4d'%('Valence electrons',self.valence_electrons[atom.element]))
         atom.steric_number += self.valence_electrons[atom.element]
-        
+
         self.display('%65s: %4d'%('Number of bonds',len(atom.bonded_atoms)))
         atom.steric_number += len(atom.bonded_atoms)
 
@@ -216,7 +219,7 @@ class Protonate:
 
         self.display('%65s: %4.1f'%('Charge(-)',atom.charge))
         atom.steric_number -= atom.charge
-        
+
         atom.steric_number = math.floor(atom.steric_number/2.0)
 
         atom.number_of_lone_pairs = atom.steric_number - len(atom.bonded_atoms) - atom.number_of_protons_to_add
@@ -240,9 +243,9 @@ class Protonate:
 
         return
 
-    
+
     def trigonal(self, atom):
-        self.display('TRIGONAL - %d bonded atoms'%(len(atom.bonded_atoms))) 
+        self.display('TRIGONAL - %d bonded atoms'%(len(atom.bonded_atoms)))
         rot_angle = math.radians(120.0)
 
         c = vector(atom1 = atom)
@@ -250,13 +253,13 @@ class Protonate:
         # 0 bonds
         if len(atom.bonded_atoms) == 0:
             pass
-            
+
         # 1 bond
         if len(atom.bonded_atoms) == 1 and atom.number_of_protons_to_add > 0:
             # Add another atom with the right angle to the first one
             a = vector(atom1 = atom, atom2 = atom.bonded_atoms[0])
             # use plane of bonded trigonal atom - e.g. arg
-            
+
             self.set_steric_number_and_lone_pairs(atom.bonded_atoms[0])
             if atom.bonded_atoms[0].steric_number == 3 and len(atom.bonded_atoms[0].bonded_atoms)>1:
                 # use other atoms bonded to the neighbour to establish the plane, if possible
@@ -265,21 +268,21 @@ class Protonate:
                     if atom.bonded_atoms[0].bonded_atoms[i] != atom:
                         other_atom_indices.append(i)
 
-                
+
                 v1 = vector(atom1 = atom, atom2 = atom.bonded_atoms[0])
-                v2 = vector(atom1 = atom.bonded_atoms[0], 
+                v2 = vector(atom1 = atom.bonded_atoms[0],
                             atom2 = atom.bonded_atoms[0].bonded_atoms[other_atom_indices[0]])
-              
+
                 axis = v1**v2
- 
+
                  # this is a trick to make sure that the order of atoms doesn't influence
                  # the final postions of added protons
                 if len(other_atom_indices)>1:
-                    v3 = vector(atom1 = atom.bonded_atoms[0], 
+                    v3 = vector(atom1 = atom.bonded_atoms[0],
                                 atom2 = atom.bonded_atoms[0].bonded_atoms[other_atom_indices[1]])
 
                     axis2 = v1**v3
-                    
+
                     if axis * axis2>0:
                         axis = axis+axis2
                     else:
@@ -294,12 +297,12 @@ class Protonate:
 
         # 2 bonds
         if len(atom.bonded_atoms) == 2 and atom.number_of_protons_to_add > 0:
-            # Add another atom with the right angle to the first two 
+            # Add another atom with the right angle to the first two
             a1 = vector(atom1 = atom, atom2 = atom.bonded_atoms[0]).rescale(1.0)
             a2 = vector(atom1 = atom, atom2 = atom.bonded_atoms[1]).rescale(1.0)
-        
+
-            new_a = -a1 - a2 
+            new_a = -a1 - a2
-            new_a = self.set_bond_distance(new_a, atom.element)            
+            new_a = self.set_bond_distance(new_a, atom.element)
             self.add_proton(atom, c+new_a)
 
 
@@ -307,20 +310,20 @@ class Protonate:
 
 
     def tetrahedral(self, atom):
-        self.display('TETRAHEDRAL - %d bonded atoms'%(len(atom.bonded_atoms))) 
+        self.display('TETRAHEDRAL - %d bonded atoms'%(len(atom.bonded_atoms)))
         rot_angle = math.radians(109.5)
 
         # sanity check
         # if atom.number_of_protons_to_add + len(atom.bonded_atoms) != 4:
-        # self.display 'Error: Attempting tetrahedral structure with %d bonds'%(atom.number_of_protons_to_add + 
+        # self.display 'Error: Attempting tetrahedral structure with %d bonds'%(atom.number_of_protons_to_add +
         #                                                                len(atom.bonded_atoms))
-        
+
         c = vector(atom1 = atom)
 
         # 0 bonds
         if len(atom.bonded_atoms) == 0:
             pass
-                  
+
         # 1 bond
         if len(atom.bonded_atoms) == 1 and atom.number_of_protons_to_add > 0:
             # Add another atom with the right angle to the first one
@@ -332,7 +335,7 @@ class Protonate:
 
         # 2 bonds
         if len(atom.bonded_atoms) == 2 and atom.number_of_protons_to_add > 0:
-            # Add another atom with the right angle to the first two             
+            # Add another atom with the right angle to the first two
             a1 = vector(atom1 = atom, atom2 = atom.bonded_atoms[0]).rescale(1.0)
             a2 = vector(atom1 = atom, atom2 = atom.bonded_atoms[1]).rescale(1.0)
 
@@ -351,16 +354,16 @@ class Protonate:
             new_a =  -a1-a2-a3
             new_a = self.set_bond_distance(new_a, atom.element)
             self.add_proton(atom, c+new_a)
-      
+
         return
 
 
     def add_proton(self, atom, position):
         # Create the new proton
-        new_H = Source.atom.Atom()
+        new_H = propka.atom.Atom()
-        new_H.setProperty(numb    = None, 
+        new_H.setProperty(numb    = None,
-                          name    = 'H%s'%atom.name[1:], 
+                          name    = 'H%s'%atom.name[1:],
-                          resName = atom.resName, 
+                          resName = atom.resName,
                           chainID = atom.chainID,
                           resNumb = atom.resNumb,
                           x       = round(position.x,3), # round of to three digimal points
@@ -381,19 +384,19 @@ class Protonate:
 
         atom.bonded_atoms.append(new_H)
         atom.number_of_protons_to_add -=1
-        atom.conformation_container.add_atom(new_H)        
+        atom.conformation_container.add_atom(new_H)
 
         # update names of all protons on this atom
-        new_H.residue_label = "%-3s%4d%2s" % (new_H.name,new_H.resNumb, new_H.chainID)        
+        new_H.residue_label = "%-3s%4d%2s" % (new_H.name,new_H.resNumb, new_H.chainID)
         no_protons = atom.count_bonded_elements('H')
         if no_protons > 1:
             i = 1
             for proton in atom.get_bonded_elements('H'):
                 proton.name = 'H%s%d'%(atom.name[1:],i)
-                proton.residue_label = "%-3s%4d%2s" % (proton.name,proton.resNumb, proton.chainID)        
+                proton.residue_label = "%-3s%4d%2s" % (proton.name,proton.resNumb, proton.chainID)
                 i+=1
 
-        
+
         self.display('added',new_H, 'to',atom)
         return
 
@@ -407,7 +410,7 @@ class Protonate:
         a = a.rescale(d)
 
         return a
-    
+
     def display(self,*text):
         if self.verbose:
             s = ''
@@ -429,11 +432,11 @@ if __name__ == '__main__':
         print('Error: Could not find \"%s\"'%filename)
         sys.exit(1)
 
-    
+
     p = Protonate()
     pdblist = pdb.readPDB(filename)
     my_protein = protein.Protein(pdblist,'test.pdb')
-    
+
     p.remove_all_hydrogen_atoms_from_protein(my_protein)
     my_protein.writePDB('before_protonation.pdb')
 

propka/run.py

@@ -0,0 +1,15 @@
+# entry point for propka script
+
+import propka.lib, propka.molecular_container
+
+def main():
+    """
+    Reads in structure files, calculates pKa values, and prints pKa files
+    """
+    # loading options, flaggs and arguments
+    options, pdbfiles = propka.lib.loadOptions()
+
+    for pdbfile in pdbfiles:
+        my_molecule = propka.molecular_container.Molecular_container(pdbfile, options)
+        my_molecule.calculate_pka()
+        my_molecule.write_pka()

propka/vector_algebra.py

@@ -1,3 +1,5 @@
+from __future__ import division
+from __future__ import print_function
 import math
 
 class vector:

propka/version.py

@@ -1,8 +1,11 @@
+from __future__ import division
+from __future__ import print_function
 import math
-import Source.lib as lib
 import sys, os
-import Source.calculations as calculations
+
-import Source.parameters
+import propka.lib as lib
+import propka.calculations as calculations
+import propka.parameters
 
 
 class version:
@@ -21,7 +24,7 @@ class version:
     def hydrogen_bond_interaction(self, group1, group2):
         return self.hydrogen_bond_interaction_model(group1, group2, self)
 
-    def calculateSideChainEnergy(self, distance, dpka_max, cutoff, weight, f_angle):    
+    def calculateSideChainEnergy(self, distance, dpka_max, cutoff, weight, f_angle):
         return self.sidechain_interaction_model(distance, dpka_max, cutoff, f_angle) # weight is ignored in 3.0 Sep07
 
     # coulomb
@@ -47,38 +50,38 @@ class version:
 
     def setup_bonding(self, molecular_container):
         return self.prepare_bonds(self.parameters, molecular_container)
- 
+
 
 
 class version_A(version):
     def __init__(self, parameters):
-        # set the calculation rutines used in this version        
+        # set the calculation rutines used in this version
         version.__init__(self, parameters)
 
         # atom naming, bonding, and protonation
-        self.molecular_preparation_method = Source.calculations.setup_bonding_and_protonation
+        self.molecular_preparation_method = propka.calculations.setup_bonding_and_protonation
-        self.prepare_bonds = Source.calculations.setup_bonding
+        self.prepare_bonds = propka.calculations.setup_bonding
 
 
         # desolvation related methods
         self.desolvation_model = calculations.radial_volume_desolvation
         self.weight_pair_method = calculations.calculatePairWeight
-        
+
         # side chain methods
-        self.sidechain_interaction_model = Source.calculations.HydrogenBondEnergy
+        self.sidechain_interaction_model = propka.calculations.HydrogenBondEnergy
-        self.hydrogen_bond_interaction_model = Source.calculations.hydrogen_bond_interaction
+        self.hydrogen_bond_interaction_model = propka.calculations.hydrogen_bond_interaction
 
         # colomb methods
-        self.electrostatic_interaction_model = Source.calculations.electrostatic_interaction
+        self.electrostatic_interaction_model = propka.calculations.electrostatic_interaction
-        self.check_coulomb_pair_method = Source.calculations.checkCoulombPair
+        self.check_coulomb_pair_method = propka.calculations.checkCoulombPair
-        self.coulomb_interaction_model = Source.calculations.CoulombEnergy
+        self.coulomb_interaction_model = propka.calculations.CoulombEnergy
 
         #backbone
-        self.backbone_interaction_model = Source.calculations.HydrogenBondEnergy
+        self.backbone_interaction_model = propka.calculations.HydrogenBondEnergy
-        self.backbone_reorganisation_method = Source.calculations.BackBoneReorganization
+        self.backbone_reorganisation_method = propka.calculations.BackBoneReorganization
 
         # exception methods
-        self.exception_check_method = Source.calculations.checkExceptions
+        self.exception_check_method = propka.calculations.checkExceptions
         return
 
     def get_hydrogen_bond_parameters(self, atom1, atom2):
@@ -91,7 +94,7 @@ class version_A(version):
             if atom.group_type in self.parameters.backbone_CO_hydrogen_bond.keys():
                 [v,c1,c2] = self.parameters.backbone_CO_hydrogen_bond[atom.group_type]
                 return [v,[c1,c2]]
-            
+
         if backbone_atom.group_type == 'BBN':
             if atom.group_type in self.parameters.backbone_NH_hydrogen_bond.keys():
                 [v,c1,c2] = self.parameters.backbone_NH_hydrogen_bond[atom.group_type]
@@ -104,7 +107,7 @@ class version_A(version):
 
 class simple_hb(version_A):
     def __init__(self, parameters):
-        # set the calculation rutines used in this version        
+        # set the calculation rutines used in this version
         version_A.__init__(self, parameters)
         print('Using simple hb model')
         return
@@ -116,23 +119,23 @@ class simple_hb(version_A):
     def get_backbone_hydrogen_bond_parameters(self, backbone_atom, atom):
         return self.parameters.hydrogen_bonds.get_value(backbone_atom.element, atom.element)
 
-    
+
 
 
 class element_based_ligand_interactions(version_A):
     def __init__(self, parameters):
-        # set the calculation rutines used in this version        
+        # set the calculation rutines used in this version
         version_A.__init__(self, parameters)
         print('Using detailed SC model!')
         return
 
     def get_hydrogen_bond_parameters(self, atom1, atom2):
-        if not 'hetatm' in [atom1.type, atom2.type]: 
+        if not 'hetatm' in [atom1.type, atom2.type]:
             # this is a protein-protein interaction
             dpka_max = self.parameters.sidechain_interaction.get_value(atom1.group_type, atom2.group_type)
             cutoff   = self.parameters.sidechain_cutoffs.get_value(atom1.group_type, atom2.group_type)
             return [dpka_max, cutoff]
-        
+
         # at least one ligand atom is involved in this interaction
         # make sure that we are using the heavy atoms for finding paramters
         elements = []
@@ -140,17 +143,17 @@ class element_based_ligand_interactions(version_A):
             if a.element == 'H': elements.append(a.bonded_atoms[0].element)
             else: elements.append(a.element)
 
-        return self.parameters.hydrogen_bonds.get_value(elements[0], elements[1])            
+        return self.parameters.hydrogen_bonds.get_value(elements[0], elements[1])
+
 
-        
     def get_backbone_hydrogen_bond_parameters(self, backbone_atom, atom):
-        if atom.type == 'atom': 
+        if atom.type == 'atom':
             # this is a backbone-protein interaction
             if backbone_atom.group_type == 'BBC' and\
                     atom.group_type in self.parameters.backbone_CO_hydrogen_bond.keys():
                 [v,c1,c2] = self.parameters.backbone_CO_hydrogen_bond[atom.group_type]
                 return [v,[c1,c2]]
-            
+
             if backbone_atom.group_type == 'BBN' and\
                     atom.group_type in self.parameters.backbone_NH_hydrogen_bond.keys():
                 [v,c1,c2] = self.parameters.backbone_NH_hydrogen_bond[atom.group_type]
@@ -168,7 +171,7 @@ class element_based_ligand_interactions(version_A):
                 print('Could not determine backbone interaction parameters for:',
                       backbone_atom,atom)
 
-            return 
+            return
 
         return None
 
@@ -176,28 +179,28 @@ class element_based_ligand_interactions(version_A):
 
 class propka30(version):
     def __init__(self, parameters):
-        # set the calculation rutines used in this version        
+        # set the calculation rutines used in this version
         version.__init__(self, parameters)
 
         # atom naming, bonding, and protonation
-        self.molecular_preparation_method = Source.calculations.setup_bonding_and_protonation_30_style
+        self.molecular_preparation_method = propka.calculations.setup_bonding_and_protonation_30_style
 
         # desolvation related methods
         self.desolvation_model = calculations.radial_volume_desolvation
         self.weight_pair_method = calculations.calculatePairWeight
-        
+
         # side chain methods
-        self.sidechain_interaction_model = Source.calculations.HydrogenBondEnergy
+        self.sidechain_interaction_model = propka.calculations.HydrogenBondEnergy
 
         # colomb methods
-        self.check_coulomb_pair_method = Source.calculations.checkCoulombPair
+        self.check_coulomb_pair_method = propka.calculations.checkCoulombPair
-        self.coulomb_interaction_model = Source.calculations.CoulombEnergy
+        self.coulomb_interaction_model = propka.calculations.CoulombEnergy
 
         #backbone
-        self.backbone_reorganisation_method = Source.calculations.BackBoneReorganization
+        self.backbone_reorganisation_method = propka.calculations.BackBoneReorganization
 
         # exception methods
-        self.exception_check_method = Source.calculations.checkExceptions
+        self.exception_check_method = propka.calculations.checkExceptions
 
 
         return

scripts/propka31.py

@@ -0,0 +1,31 @@
+#!/usr/bin/env python
+
+# This is the original propka script. However, this distribute-based
+# installation moved the main() function into propka.run.main and just
+# generates a script called propka31 from the setup.py installation
+# script. You should not need to use this script.
+#
+# (Also note that there can be import problems because the script name
+# is the same as the module name; that's why the new script is called
+# propka31.)
+
+import propka.lib, propka.molecular_container
+
+def main():
+    """
+    Reads in structure files, calculates pKa values, and prints pKa files
+    """
+    # loading options, flaggs and arguments
+    options, pdbfiles = propka.lib.loadOptions()
+
+    for pdbfile in pdbfiles:
+        my_molecule = propka.molecular_container.Molecular_container(pdbfile, options)
+        my_molecule.calculate_pka()
+        my_molecule.write_pka()
+
+
+if __name__ == '__main__':
+    #import cProfile
+    #cProfile.run('main()',sort=1)
+    main()
+

setup.py

@@ -0,0 +1,51 @@
+# PROPKA 3.1
+#
+#
+# setuptools installation of  PROPKA 3.1
+
+import ez_setup
+ez_setup.use_setuptools()
+from setuptools import setup, find_packages
+
+VERSION = "3.1"
+
+setup(name="PROPKA",
+      version=VERSION,
+      description="Heuristic pKa calculations with ligands",
+      long_description="""
+PROPKA predicts the pKa values of ionizable groups in proteins (version 3.0) and
+protein-ligand complexes (version 3.1) based on the 3D structure.
+
+For proteins without ligands both version should produce the same result.
+
+The method is described in the following papers, which you should cite
+in publications:
+
+* Sondergaard, Chresten R., Mats HM Olsson, Michal Rostkowski, and Jan
+  H. Jensen. "Improved Treatment of Ligands and Coupling Effects in
+  Empirical Calculation and Rationalization of pKa Values." Journal of
+  Chemical Theory and Computation 7, no. 7 (2011): 2284-2295.
+
+* Olsson, Mats HM, Chresten R. Sondergaard, Michal Rostkowski, and Jan
+  H. Jensen. "PROPKA3: consistent treatment of internal and surface
+  residues in empirical pKa predictions." Journal of Chemical Theory
+  and Computation 7, no. 2 (2011): 525-537.
+
+See http://propka.ki.ku.dk/ for the PROPKA web server,
+using the tutorial http://propka.ki.ku.dk/~luca/wiki/index.php/PROPKA_3.1_Tutorial .
+""",
+      author="Jan H. Jensen",
+      author_email="jhjensen@chem.ku.dk",
+      license="",
+      url="http://propka.ki.ku.dk/",
+      keywords="science",
+      packages=find_packages(exclude=['scripts']),
+      package_data = {'propka': ['*.dat', '*.cfg']},
+      #scripts = ["scripts/propka31.py"],  # use entry point below
+      entry_points = {
+        'console_scripts': [
+            'propka31 = propka.run:main',
+            ],
+        },
+      zip_safe=True,
+)