Package: pdb2pqr | Debian Sources

Package: pdb2pqr / 2.1.1+dfsg-7+deb11u1

python3_sundry.patch
Description: Use 2to3 to port from Python2 to Python3
Bug-Debian: https://bugs.debian.org/937262
Author: Andreas Tille <tille@debian.org>
Last-Update: Tue, 15 Feb 2022 07:03:38 -0500

--- a/main_cgi.py
+++ b/main_cgi.py
@@ -53,7 +53,7 @@
 import re
 import sys
 import string
-from StringIO import StringIO
+from io import StringIO
 #import tempfile
 #from src import pdb
 
@@ -94,16 +94,16 @@
     """
     if jobid:
         if have_opal:
-            print "Location: querystatus.cgi?jobid=%s&typeofjob=opal\n" % (jobid)
+            print("Location: querystatus.cgi?jobid=%s&typeofjob=opal\n" % (jobid))
         else:
-            print "Location: querystatus.cgi?jobid=%s&typeofjob=local\n" % (jobid)
+            print("Location: querystatus.cgi?jobid=%s&typeofjob=local\n" % (jobid))
 
     #print "Content-type: text/html\n"
-    print "<HTML>"
-    print "<HEAD>"
-    print "\t<TITLE>%s</TITLE>" % pagetitle
-    print "\t<link rel=\"stylesheet\" href=\"%s\" type=\"text/css\">\n" % STYLESHEET
-    print "</HEAD>"
+    print("<HTML>")
+    print("<HEAD>")
+    print("\t<TITLE>%s</TITLE>" % pagetitle)
+    print("\t<link rel=\"stylesheet\" href=\"%s\" type=\"text/css\">\n" % STYLESHEET)
+    print("</HEAD>")
     return
 
 def redirector(name, weboptions):
@@ -135,7 +135,7 @@
     if 'selectedExtensions' in analiticsDict:
         analiticsDict['selectedExtensions'] = ' '.join(analiticsDict['selectedExtensions'])
     
-    options = ','.join(str(k)+':'+str(v) for k,v in analiticsDict.iteritems())
+    options = ','.join(str(k)+':'+str(v) for k,v in analiticsDict.items())
     events['options']=options
 
     eventsScriptString = ''
@@ -188,31 +188,31 @@
         #These are included in has_key(), __contains__(), and __getitem__() calls.
         self.otheroptions = {}
         
-        self.runoptions['debump'] = form.has_key("DEBUMP")
-        self.runoptions['opt'] = form.has_key("OPT")
+        self.runoptions['debump'] = "DEBUMP" in form
+        self.runoptions['opt'] = "OPT" in form
         
-        if form.has_key('FF'):
+        if 'FF' in form:
             self.ff = form["FF"].value.lower()
         else:
             raise WebOptionsError('Force field type missing from form.')
         
-        if form.has_key("PDBID") and form["PDBID"].value and form["PDBSOURCE"].value == 'ID':
+        if "PDBID" in form and form["PDBID"].value and form["PDBSOURCE"].value == 'ID':
             self.pdbfile = utilities.getPDBFile(form["PDBID"].value)
             if self.pdbfile is None:
                 raise WebOptionsError('The pdb ID provided is invalid.')
             self.pdbfilestring = self.pdbfile.read()
             self.pdbfile = StringIO(self.pdbfilestring)
             self.pdbfilename = form["PDBID"].value
-        elif form.has_key("PDB") and form["PDB"].filename and form["PDBSOURCE"].value == 'UPLOAD':
+        elif "PDB" in form and form["PDB"].filename and form["PDBSOURCE"].value == 'UPLOAD':
             self.pdbfilestring = form["PDB"].value
             self.pdbfile = StringIO(self.pdbfilestring)
             self.pdbfilename = sanitizeFileName(form["PDB"].filename)
         else:
             raise WebOptionsError('You need to specify a pdb ID or upload a pdb file.')
             
-        if form.has_key("PKACALCMETHOD"):
+        if "PKACALCMETHOD" in form:
             if form["PKACALCMETHOD"].value != 'none':
-                if not form.has_key('PH'):
+                if 'PH' not in form:
                     raise WebOptionsError('Please provide a pH value.')
                 
                 phHelp = 'Please choose a pH between 0.0 and 14.0.'
@@ -237,11 +237,11 @@
                                                           'sdie': 80,
                                                           'pairene': 1.0}
                  
-        self.otheroptions['apbs'] = form.has_key("INPUT")
-        self.otheroptions['whitespace'] = form.has_key("WHITESPACE")
+        self.otheroptions['apbs'] = "INPUT" in form
+        self.otheroptions['whitespace'] = "WHITESPACE" in form
         
         if self.ff == 'user':
-            if form.has_key("USERFF") and form["USERFF"].filename:
+            if "USERFF" in form and form["USERFF"].filename:
                 self.userfffilename = sanitizeFileName(form["USERFF"].filename)
                 self.userffstring = form["USERFF"].value
                 self.runoptions['userff'] = StringIO(form["USERFF"].value)
@@ -249,7 +249,7 @@
                 text = "A force field file must be provided if using a user created force field."
                 raise WebOptionsError(text)
                 
-            if form.has_key("USERNAMES") and form["USERNAMES"].filename:
+            if "USERNAMES" in form and form["USERNAMES"].filename:
                 self.usernamesfilename = sanitizeFileName(form["USERNAMES"].filename)
                 self.usernamesstring = form["USERNAMES"].value
                 self.runoptions['usernames'] = StringIO(form["USERNAMES"].value)
@@ -257,20 +257,20 @@
                 text = "A names file must be provided if using a user created force field."
                 raise WebOptionsError(text)
             
-        if form.has_key("FFOUT") and form["FFOUT"].value != "internal":
+        if "FFOUT" in form and form["FFOUT"].value != "internal":
             self.runoptions['ffout'] = form["FFOUT"].value
                 
-        self.runoptions['chain'] = form.has_key("CHAIN")
-        self.runoptions['typemap'] = form.has_key("TYPEMAP")
-        self.runoptions['neutraln'] = form.has_key("NEUTRALN")
-        self.runoptions['neutralc'] = form.has_key("NEUTRALC")
-        self.runoptions['drop_water'] = form.has_key("DROPWATER")
+        self.runoptions['chain'] = "CHAIN" in form
+        self.runoptions['typemap'] = "TYPEMAP" in form
+        self.runoptions['neutraln'] = "NEUTRALN" in form
+        self.runoptions['neutralc'] = "NEUTRALC" in form
+        self.runoptions['drop_water'] = "DROPWATER" in form
         
         if (self.runoptions['neutraln'] or self.runoptions['neutraln']) and \
             self.ff != 'parse':
             raise WebOptionsError('Neutral N-terminus and C-terminus require the PARSE forcefield.')
         
-        if form.has_key("LIGAND") and form['LIGAND'].filename:
+        if "LIGAND" in form and form['LIGAND'].filename:
             self.ligandfilename=sanitizeFileName(form["LIGAND"].filename)
             ligandfilestring = form["LIGAND"].value
             # for Windows and Mac style newline compatibility for pdb2pka
@@ -313,16 +313,16 @@
         options['pdb'] = self.pdbfilename
         
         #propkaOptions is redundant.
-        if options.has_key('ph_calc_options'):
+        if 'ph_calc_options' in options:
             del options['ph_calc_options']
         
-        if options.has_key('ligand'):
+        if 'ligand' in options:
             options['ligand'] = self.ligandfilename
             
-        if options.has_key('userff'):
+        if 'userff' in options:
             options['userff'] = self.userfffilename
             
-        if options.has_key('usernames'):
+        if 'usernames' in options:
             options['usernames'] = self.usernamesfilename
         
         return options
@@ -437,24 +437,24 @@
 
     try:
         resp=appServicePort.launchJob(req)
-    except Exception, e:
+    except Exception as e:
         printHeader("PDB2PQR Job Submission - Error")
-        print "<BODY>\n<P>"
-        print "There was an error with your job submission<br>"
-        print "</P>"
-        print "<script type=\"text/javascript\">"
-        print "var gaJsHost = ((\"https:\" == document.location.protocol) ? \"https://ssl.\" : \"http://www.\");"
-        print "document.write(unescape(\"%3Cscript src=\'\" + gaJsHost + \"google-analytics.com/ga.js\' type=\'text/javascript\'%3E%3C/script%3E\"));"
-        print "</script>"
-        print "<script type=\"text/javascript\">"
-        print "try {"
-        print "var pageTracker = _gat._getTracker(\"UA-11026338-3\");"
+        print("<BODY>\n<P>")
+        print("There was an error with your job submission<br>")
+        print("</P>")
+        print("<script type=\"text/javascript\">")
+        print("var gaJsHost = ((\"https:\" == document.location.protocol) ? \"https://ssl.\" : \"http://www.\");")
+        print("document.write(unescape(\"%3Cscript src=\'\" + gaJsHost + \"google-analytics.com/ga.js\' type=\'text/javascript\'%3E%3C/script%3E\"));")
+        print("</script>")
+        print("<script type=\"text/javascript\">")
+        print("try {")
+        print("var pageTracker = _gat._getTracker(\"UA-11026338-3\");")
         for key in weboptions:
-            print "pageTracker._trackPageview(\"/main_cgi/has_%s_%s.html\");" % (key, weboptions[key])
-        print "pageTracker._trackPageview();"
-        print "} catch(err) {}</script>"
-        print "</BODY>"
-        print "</HTML>"
+            print("pageTracker._trackPageview(\"/main_cgi/has_%s_%s.html\");" % (key, weboptions[key]))
+        print("pageTracker._trackPageview();")
+        print("} catch(err) {}</script>")
+        print("</BODY>")
+        print("</HTML>")
         sys.exit(2)
     
     try:
@@ -475,7 +475,7 @@
         pdb2pqrOpalJobIDFile.write(resp._jobID)
         pdb2pqrOpalJobIDFile.close()
         
-        print redirector(name, weboptions)
+        print(redirector(name, weboptions))
         
         # Recording CGI run information for PDB2PQR Opal
         pdb2pqrOpalLogFile = open('%s%s%s/pdb2pqr_log' % (INSTALLDIR, TMPDIR, name), 'w')
@@ -483,8 +483,8 @@
                                  str(os.environ["REMOTE_ADDR"]))
         pdb2pqrOpalLogFile.close()
 
-    except StandardError, details:
-        print details
+    except Exception as details:
+        print(details)
         createError(name, details)
 
 def handleNonOpal(weboptions):
@@ -500,7 +500,7 @@
         text = "Unable to find PDB file - Please make sure this is "
         text += "a valid PDB file ID!"
         #print "Content-type: text/html\n"
-        print text
+        print(text)
         sys.exit(2)
     elif dummyprot.numAtoms() > MAXATOMS and weboptions["opt"] == True:
         text = "<HTML><HEAD>"
@@ -527,7 +527,7 @@
         text += "} catch(err) {}</script>"
         text += "</BODY></HTML>"
         #print "Content-type: text/html\n"
-        print text
+        print(text)
         sys.exit(2)
 
     try:
@@ -558,7 +558,7 @@
 
         pid = os.fork()
         if pid:
-            print redirector(name, weboptions)
+            print(redirector(name, weboptions))
             sys.exit()
         else:
             currentdir = os.getcwd()
@@ -646,10 +646,10 @@
 
     #TODO: Better error reporting.
     #Also, get forked job to properly write error status on failure.
-    except StandardError, details:
+    except Exception as details:
     #except StandardError as details:
-        print traceback.format_exc()
-        print sys.exc_info()[0]
+        print(traceback.format_exc())
+        print(sys.exc_info()[0])
         #print details
         createError(name, details)
 
@@ -657,7 +657,7 @@
     """
         Main driver for running PDB2PQR from a web page
     """
-    print "Content-type: text/html\n"
+    print("Content-type: text/html\n")
     import cgi
     import cgitb
 
--- a/pdb2pka/NEWligand_topology.py
+++ b/pdb2pka/NEWligand_topology.py
@@ -7,13 +7,13 @@
 import numpy
     
 #from sets import Set
-from ligandclean.trial_templates import *
-from ligandclean.lookuptable import *
+from .ligandclean.trial_templates import *
+from .ligandclean.lookuptable import *
 try:
-    from substruct import Algorithms
+    from .substruct import Algorithms
 except ImportError:
     txt = "Cannot import Algorithms, this may be the result of disabling pdb2pka at configure stage!"    
-    raise ImportError, txt
+    raise ImportError(txt)
 from types import *
 
 def length(vector):
@@ -51,7 +51,7 @@
             # Get the likely types from names
             #
             trivial_types=['N','O','C','H']
-            for atom in self.atoms.keys():
+            for atom in list(self.atoms.keys()):
                 if atom[0] in trivial_types:
                     if atom[0]!='H': # Get rid of all the hydrogens
                         self.atoms[atom]['type']=atom[0]
@@ -61,9 +61,9 @@
             # First approximation: Anything closer than 2.0 A is bonded
             #
             self.dists={}
-            for atom1 in self.atoms.keys():
+            for atom1 in list(self.atoms.keys()):
                 self.dists[atom1]={}
-                for atom2 in self.atoms.keys():
+                for atom2 in list(self.atoms.keys()):
                     if atom1==atom2:
                         continue
                     #
@@ -79,7 +79,7 @@
             #
             # Get the torsion angles
             #
-            atoms=self.atoms.keys()
+            atoms=list(self.atoms.keys())
             atoms.sort()
             for atom in atoms:
                 self.atoms[atom]['torsions']=self.get_torsions(atom)
@@ -128,7 +128,7 @@
             #
             # Get the torsion angles
             #
-            atoms=self.atoms.keys()
+            atoms=list(self.atoms.keys())
             atoms.sort()
             for atom in atoms:
                 self.atoms[atom]['torsions']=self.get_torsions(atom)
@@ -149,7 +149,7 @@
         # + determine their likely order (e.g. single, double, or triple)
         #
         ambs={}
-        atoms=self.atoms.keys()
+        atoms=list(self.atoms.keys())
         for atom_name in atoms:
             bonds=self.atoms[atom_name]['bonds']
             atype=self.atoms[atom_name]['type']
@@ -172,13 +172,13 @@
         #
         valences={'C':4,'O':2,'N':3}
         for atom in self.atoms:
-            print atom, self.atoms[atom]
+            print(atom, self.atoms[atom])
         #
         # ok, now it gets hairy
         #
-        print
-        print 'Guessing sybyl atom types'
-        for atom in self.atoms.keys():
+        print()
+        print('Guessing sybyl atom types')
+        for atom in list(self.atoms.keys()):
             stype=None
             at=self.atoms[atom]
             #
@@ -218,7 +218,7 @@
         # Do some postchecks
         # - right now only for Carboxylic acids
         #
-        for atom in self.atoms.keys():
+        for atom in list(self.atoms.keys()):
             at=self.atoms[atom]
             if at['sybylType']=='C.2':
                 #
@@ -239,11 +239,11 @@
         #
         # All Done
         #
-        atoms=self.atoms.keys()
+        atoms=list(self.atoms.keys())
         atoms.sort()
-        print '\nFinal Sybyl type results'
+        print('\nFinal Sybyl type results')
         for atom in atoms:
-            print atom,self.atoms[atom]['sybylType']
+            print(atom,self.atoms[atom]['sybylType'])
         return
 
     #
@@ -273,7 +273,7 @@
         tps.sort() # To agree with dictionary layout
         best_fit=2.00
         best_type=None
-        for bond in bond_props.keys():
+        for bond in list(bond_props.keys()):
             if bond[0]==tps[0] and bond[-1]==tps[1]:
                 if abs(dist-bond_props[bond][0])<best_fit:
                     best_fit=abs(dist-bond_props[bond][0])
@@ -322,7 +322,7 @@
         # Make the lines for the pdb2pqr definition
         #
         self.numbers={}
-        atoms=self.atoms.keys()
+        atoms=list(self.atoms.keys())
         atoms.sort()
         number=0
         for atom in atoms:
@@ -429,12 +429,12 @@
         """#
         # Look for simple substructures that would be titratable groups in the ligand
         #"""
-        atoms=self.atoms.keys()
+        atoms=list(self.atoms.keys())
         #
         # ring detection (including deleting redundancies & sorting issues)
         ring_list = []
         tmp=[]
-        for atom in self.atoms.keys():
+        for atom in list(self.atoms.keys()):
             temp_ring_list = []
             tmp.append(self.ring_detection(atom))
         #
@@ -462,7 +462,7 @@
         #
         #
         # assigning ring attribute for every ring atom
-        for at in self.atoms.keys():
+        for at in list(self.atoms.keys()):
             self.atoms[at]['in_ring'] = 0
         for rring in sorted_ring_list:
             for current_atom in rring:
@@ -556,7 +556,7 @@
         dictCounter = 0
         dict_of_matched_lig_fragments = {}
         matched_lig_template = []
-        for current_template in templates.keys():
+        for current_template in list(templates.keys()):
             #print "MATCHING", current_template
             #print "######## ########"
             output = match(templates[current_template],atoms)
@@ -571,7 +571,7 @@
                 temp_templateList = list(set(templateList))
                 matchedLigAtoms =[]
                 temptemp = []
-                if len(temp_templateList) == len(ligList) and (len(ee) == len(templates[current_template].keys())):
+                if len(temp_templateList) == len(ligList) and (len(ee) == len(list(templates[current_template].keys()))):
                     for xxee in ee:
                         # output[1][xxee] => maching pair of ligand and template atoms
                         # temporary depostion
@@ -631,10 +631,10 @@
                                         # loop over all entries
                                         for possiblyredundantentries in NonRedundantCliques:
                                             if set(possiblyredundantentries).issubset(set(xxxx)):
-                                                print NonRedundantCliques
+                                                print(NonRedundantCliques)
                                                 NonRedundantCliques.remove(possiblyredundantentries)
                                                 NonRedundantCliques.append(xxxx)
-                                                print NonRedundantCliques
+                                                print(NonRedundantCliques)
                                             elif set(xxxx).issubset(set(possiblyredundantentries)):
                                                 #print "found subset which is not added to the list"
                                                 pass
@@ -665,17 +665,17 @@
          
         for allCl in dict_of_matched_lig_fragments:
             if dict_of_matched_lig_fragments[allCl]['matchedligatoms'] == NonRedundantCliques[0]:
-                print "WE MATCHED", dict_of_matched_lig_fragments[allCl]['templatename']
-                print "matchedligatoms            : ", dict_of_matched_lig_fragments[allCl]['matchedligatoms']
-                print "type                       : ", dict_of_matched_lig_fragments[allCl]['type']
-                print "modelpka                   : ", dict_of_matched_lig_fragments[allCl]['modelpka']
-                print "titratableatoms            : ", dict_of_matched_lig_fragments[allCl]['titratableatoms']
-                print "matching atoms             : ", dict_of_matched_lig_fragments[allCl]['matching_atoms']
+                print("WE MATCHED", dict_of_matched_lig_fragments[allCl]['templatename'])
+                print("matchedligatoms            : ", dict_of_matched_lig_fragments[allCl]['matchedligatoms'])
+                print("type                       : ", dict_of_matched_lig_fragments[allCl]['type'])
+                print("modelpka                   : ", dict_of_matched_lig_fragments[allCl]['modelpka'])
+                print("titratableatoms            : ", dict_of_matched_lig_fragments[allCl]['titratableatoms'])
+                print("matching atoms             : ", dict_of_matched_lig_fragments[allCl]['matching_atoms'])
         # re-run matching to get mutiple titratable sites?
         
         # TJD: This is to resolve the bug fix when allCl is None
         if dict_of_matched_lig_fragments != {}:
-            print dict_of_matched_lig_fragments[allCl]
+            print(dict_of_matched_lig_fragments[allCl])
             return dict_of_matched_lig_fragments[allCl]
         else:
             return {}
--- a/pdb2pka/apbs.py
+++ b/pdb2pka/apbs.py
@@ -20,11 +20,11 @@
     #
     # We need _apbslib.so and apbslib.py
     #
-    print
-    print 'Missing libraries for interfacing with APBS'
-    print
-    print 'You need to build APBS with Python support using the CMake variable -DENABLE_PYTHON=ON.'
-    print
+    print()
+    print('Missing libraries for interfacing with APBS')
+    print()
+    print('You need to build APBS with Python support using the CMake variable -DENABLE_PYTHON=ON.')
+    print()
     sys.exit(0)
 
 Python_kb = 1.3806581e-23
@@ -53,7 +53,7 @@
         """
             Return the error message
         """
-        return `self.value`
+        return repr(self.value)
 
 class runAPBS:
 
@@ -226,7 +226,7 @@
     
             routines.write("Preparing to run %d PBE calculations. \n" % self.nosh.ncalc)
     
-            for icalc in xrange(self.nosh.ncalc):
+            for icalc in range(self.nosh.ncalc):
                 sys.stdout.write("---------------------------------------------\n")
                 self.calc = NOsh_getCalc(self.nosh, icalc)
                 self.mgparm = self.calc.mgparm
--- a/pdb2pka/charge_mon.py
+++ b/pdb2pka/charge_mon.py
@@ -1,4 +1,4 @@
-from Tkinter import *
+from tkinter import *
 
 class charge_mon(Frame):
 
@@ -56,13 +56,13 @@
         self.cv.create_text(0,self.calc,text=self.text,anchor='nw')
         charges={}
         for resnum,atomname,charge in charge_list:
-            if not charges.has_key(resnum):
+            if resnum not in charges:
                 charges[resnum]=[]
             charges[resnum].append(charge)
         #
         # Sum all charges
         #
-        for res in charges.keys():
+        for res in list(charges.keys()):
             non_zero=None
             sum=0.0
             for crg in charges[res]:
@@ -77,7 +77,7 @@
         #
         #
         later=[]
-        for resid in charges.keys():
+        for resid in list(charges.keys()):
             x_count=self.res_pos[resid]
             if charges[resid] is None:
                 fill='white'
@@ -98,7 +98,7 @@
         #
         for x_count,text,resid in later:
             self.cv.create_text(x_count,self.calc,text=text,anchor='nw',fill='black')
-            print '!!Wrong charge: %s %s' %(text,str(resid))
+            print('!!Wrong charge: %s %s' %(text,str(resid)))
         #
         # Update and increment row
         #
--- a/pdb2pka/example.py
+++ b/pdb2pka/example.py
@@ -46,7 +46,7 @@
         """
             Return the error message
         """
-        return `self.value`
+        return repr(self.value)
 
 def getUnitConversion():
     """
@@ -97,7 +97,7 @@
 
     if not parseInputFromString(nosh, INPUT):
         stderr.write("main:  Error while parsing input file.\n")
-        raise APBSError, "Error occurred!"
+        raise APBSError("Error occurred!")
 
     # Load the molecules using Valist_load routine
 
@@ -145,14 +145,14 @@
  
     sys.stdout.write("Preparing to run %d PBE calculations. \n" % nosh.ncalc)
    
-    for icalc in xrange(nosh.ncalc):
+    for icalc in range(nosh.ncalc):
         sys.stdout.write("---------------------------------------------\n")
         calc = NOsh_getCalc(nosh, icalc)
         mgparm = calc.mgparm
         pbeparm = calc.pbeparm
         if calc.calctype != 0:
             sys.stderr.write("main:  Only multigrid calculations supported!\n")
-            raise APBSError, "Only multigrid calculations supported!"
+            raise APBSError("Only multigrid calculations supported!")
 
         for k in range(0, nosh.nelec):
             if NOsh_elec2calc(nosh,k) >= icalc:
@@ -171,7 +171,7 @@
               alist, dielXMap, dielYMap, dielZMap, kappaMap, chargeMap, 
               pmgp, pmg) != 1:
             sys.stderr.write("Error setting up MG calculation!\n")
-            raise APBSError, "Error setting up MG calculation!"
+            raise APBSError("Error setting up MG calculation!")
 	
         # Print problem parameters 
 
@@ -183,13 +183,13 @@
 
         if solveMG(nosh, thispmg, mgparm.type) != 1:
             stderr.write("Error solving PDE! \n")
-            raise APBSError, "Error Solving PDE!"
+            raise APBSError("Error Solving PDE!")
 
         # Set partition information : Routine setPartMG
 
         if setPartMG(nosh, mgparm, thispmg) != 1:
             sys.stderr.write("Error setting partition info!\n")
-            raise APBSError, "Error setting partition info!"
+            raise APBSError("Error setting partition info!")
 	
         ret, totEnergy[icalc] = energyMG(nosh, icalc, thispmg, 0,
                                          totEnergy[icalc], 0.0, 0.0, 0.0)
@@ -212,7 +212,7 @@
     if nosh.nprint > 0:
         sys.stdout.write("---------------------------------------------\n")
         sys.stdout.write("PRINT STATEMENTS\n")
-    for iprint in xrange(nosh.nprint):
+    for iprint in range(nosh.nprint):
         if NOsh_printWhat(nosh, iprint) == NPT_ENERGY:
             printEnergy(com, nosh, totEnergy, iprint)
         elif NOsh_printWhat(nosh, iprint) == NPT_FORCE:
@@ -283,7 +283,7 @@
 
     # Print out the number of elec statements
 
-    print "Number of elecs: ", len(input.elecs)
+    print("Number of elecs: ", len(input.elecs))
 
     # Let's set the dielectric in the second elec statement
 
@@ -295,4 +295,4 @@
 
     # And print the results!
 
-    print "Now we have: ", potentials
+    print("Now we have: ", potentials)
--- a/pdb2pka/graph_cut/create_titration_output.py
+++ b/pdb2pka/graph_cut/create_titration_output.py
@@ -2,7 +2,7 @@
 
 def create_output(path, curves):
     file_name_template = "{}_{}_{}.csv"
-    for key, curve in curves.items():
+    for key, curve in list(curves.items()):
         file_name = file_name_template.format(*key)
         file_path = os.path.join(path, file_name)
 
--- a/pdb2pka/graph_cut/graph.py
+++ b/pdb2pka/graph_cut/graph.py
@@ -21,7 +21,7 @@
         self._build_nodes()
 
         #Create edges going in and out of S and T.
-        for key, v in self.pc.residue_variables.items():
+        for key, v in list(self.pc.residue_variables.items()):
             prot_instance = v.instances["PROTONATED"]
             prot_capacity = prot_instance.energyNF / 2.0
             prot_node = key+("PROTONATED",)
@@ -39,7 +39,7 @@
                 self.DG.add_edge(deprot_node, "T", capacity=deprot_capacity)
 
         #Create all interaction energy edges.
-        for p, q in combinations(iter(self.pc.residue_variables.items()),2):
+        for p, q in combinations(iter(list(self.pc.residue_variables.items())),2):
             p_key, p_residue = p
             q_key, q_residue = q
 
@@ -91,7 +91,7 @@
         """Creates a map of residues to instances based on the """
         labeling = {}
         uncertain = []
-        for key, v in self.pc.residue_variables.items():
+        for key, v in list(self.pc.residue_variables.items()):
             prot_node = key+("PROTONATED",)
             deprot_node = key+("DEPROTONATED",)
 
--- a/pdb2pka/graph_cut/protein_complex.py
+++ b/pdb2pka/graph_cut/protein_complex.py
@@ -14,7 +14,7 @@
                            "HIS": {"model_pka": 6.6, "ionizable": 1, "tautomers":{"protonated":("1", "2"), "deprotonated":("1+2",)}},
                           }
 all_tautomers = set()
-for residue_data in titratable_residue_data.values():
+for residue_data in list(titratable_residue_data.values()):
     all_tautomers.update(residue_data["tautomers"]["deprotonated"])
     all_tautomers.update(residue_data["tautomers"]["protonated"])
 state_to_tautomer = {}
@@ -77,8 +77,8 @@
         return self.instances.get(state)
     def get_prot_and_deprot_instances(self):
         #We don't want the placeholder instances
-        prot = [value for key, value in self.instances.items() if "PROTONATED" not in key and value.protonated]
-        deprot = [value for key, value in self.instances.items() if "PROTONATED" not in key and not value.protonated]
+        prot = [value for key, value in list(self.instances.items()) if "PROTONATED" not in key and value.protonated]
+        deprot = [value for key, value in list(self.instances.items()) if "PROTONATED" not in key and not value.protonated]
         return prot, deprot
     def __hash__(self):
         return hash(self.name)
@@ -180,7 +180,7 @@
         """After we have loaded the backgroind and desolvation files we need to
         update the consolidated instances to use the minimum energy of the
         instances they are meant to replace."""
-        for key, residue in self.residue_variables.items():
+        for key, residue in list(self.residue_variables.items()):
             name = key[0]
             if name == 'HIS':
                 continue
@@ -210,7 +210,7 @@
         handled_interaction_pairs = set()
 
         #See https://docs.python.org/2/library/itertools.html#itertools.combinations
-        for v, w in combinations(iter(self.residue_variables.items()),2):
+        for v, w in combinations(iter(list(self.residue_variables.items())),2):
             v_key, v_residue = v
             w_key, w_residue = w
             v_name = v_key[0]
@@ -251,7 +251,7 @@
 
         #Now handle HIS.
         #See https://docs.python.org/2/library/itertools.html#itertools.permutations
-        for v, w in permutations(iter(self.residue_variables.items()),2):
+        for v, w in permutations(iter(list(self.residue_variables.items())),2):
             his_key, his_residue = v
             other_key, other_residue = w
             his_name = his_key[0]
@@ -292,12 +292,12 @@
         #Clean out unused interaction energies
         self.drop_interaction_pairs(handled_interaction_pairs)
 
-        for key, residue in self.residue_variables.items():
+        for key, residue in list(self.residue_variables.items()):
             name = key[0]
             if name == 'HIS':
                 continue
 
-            residue.instances = OrderedDict((k,v) for k,v in residue.instances.items() if "PROTONATED" in k)
+            residue.instances = OrderedDict((k,v) for k,v in list(residue.instances.items()) if "PROTONATED" in k)
 
 
     def divide_his(self):
@@ -373,7 +373,7 @@
         #  HId <-> HIe.)
 
         #See https://docs.python.org/2/library/itertools.html#itertools.product
-        for v, w in product(iter(self.residue_variables.items()), repeat=2):
+        for v, w in product(iter(list(self.residue_variables.items())), repeat=2):
             his_key, his_residue = v
             other_key, other_residue = w
             his_name, his_chain, his_location = his_key
@@ -535,7 +535,7 @@
                     ph_multiplier += 1.0
 
             #See https://docs.python.org/2/library/itertools.html#itertools.combinations
-            for v, w in combinations(iter(self.residue_variables.items()), 2):
+            for v, w in combinations(iter(list(self.residue_variables.items())), 2):
                 v_key, v_residue = v
                 w_key, w_residue = w
                 is_same_his = (v_key[1:] == w_key[1:] and v_key[0] in ("HId", "HIe") and w_key[0] in ("HId", "HIe"))
@@ -544,12 +544,12 @@
                     if labeling[v_residue].protonated and labeling[w_residue].protonated:
                         ph_multiplier -= 2.0
 
-        for v in self.residue_variables.values():
+        for v in list(self.residue_variables.values()):
             v_instance = labeling[v]
             energy += v_instance.energyNF if normal_form else v_instance.energy
 
         ie = self.normalized_interaction_energies if normal_form else self.interaction_energies
-        for v, w in combinations(iter(self.residue_variables.values()),2):
+        for v, w in combinations(iter(list(self.residue_variables.values())),2):
             v_instance = labeling[v]
             w_instance = labeling[w]
             energy += ie[v_instance, w_instance]
@@ -562,7 +562,7 @@
         return energy
 
     def instance_interaction_energy(self, instance, labeling, interaction_energy):
-        return sum(interaction_energy.get((instance, labeling[x]),0) for x in self.residue_variables.values())
+        return sum(interaction_energy.get((instance, labeling[x]),0) for x in list(self.residue_variables.values()))
 
     def evaluate_energy_diff(self, residue, labeling, normal_form = False):
         """Returns the total energy difference between the deprontated and protonated
@@ -572,12 +572,12 @@
 
         prot_instance = residue.instances["PROTONATED"]
         prot_energy = 0
-        for x in self.residue_variables.values():
+        for x in list(self.residue_variables.values()):
             prot_energy += ie.get((prot_instance, labeling[x]),0.0)
         prot_energy += prot_instance.energyNF if normal_form else prot_instance.energy
 
         deprot_instance = residue.instances["DEPROTONATED"]
-        deprot_energy = sum(ie.get((deprot_instance, labeling[x]),0.0) for x in self.residue_variables.values())
+        deprot_energy = sum(ie.get((deprot_instance, labeling[x]),0.0) for x in list(self.residue_variables.values()))
         deprot_energy += deprot_instance.energyNF if normal_form else deprot_instance.energy
 
         return prot_energy - deprot_energy
@@ -644,8 +644,8 @@
             EHd1 = hid_prot_instance.energy
             EHd0 = hid_deprot_instance.energy
 
-        sum_ErHd_a1 = sum(ie.get((labeling_copy[x], hid_prot_instance), 0) for x in self.residue_variables.values())
-        sum_ErHd_a0 = sum(ie.get((labeling_copy[x], hid_deprot_instance), 0) for x in self.residue_variables.values())
+        sum_ErHd_a1 = sum(ie.get((labeling_copy[x], hid_prot_instance), 0) for x in list(self.residue_variables.values()))
+        sum_ErHd_a0 = sum(ie.get((labeling_copy[x], hid_deprot_instance), 0) for x in list(self.residue_variables.values()))
 
         if normal_form:
             EHe1 = hie_prot_instance.energyNF
@@ -654,8 +654,8 @@
             EHe1 = hie_prot_instance.energy
             EHe0 = hie_deprot_instance.energy
 
-        sum_ErHe_a1 = sum(ie.get((labeling_copy[x], hie_prot_instance), 0) for x in self.residue_variables.values())
-        sum_ErHe_a0 = sum(ie.get((labeling_copy[x], hie_deprot_instance), 0) for x in self.residue_variables.values())
+        sum_ErHe_a1 = sum(ie.get((labeling_copy[x], hie_prot_instance), 0) for x in list(self.residue_variables.values()))
+        sum_ErHe_a0 = sum(ie.get((labeling_copy[x], hie_deprot_instance), 0) for x in list(self.residue_variables.values()))
 
         hsp_hse = EHd1 - EHd0 + sum_ErHd_a1 - sum_ErHd_a0
         hsp_hsd = EHe1 - EHe0 + sum_ErHe_a1 - sum_ErHe_a0
@@ -670,7 +670,7 @@
         self.normalized_constant_energy = 0.0
 
         #Add the pH to the instances.
-        for residue in self.residue_variables.values():
+        for residue in list(self.residue_variables.values()):
             residue.instances["DEPROTONATED"].energyNF = residue.instances["DEPROTONATED"].energy
             # TODO - I'm not sure why this is +pH....
             residue.instances["PROTONATED"].energyNF = residue.instances["PROTONATED"].energy + math.log(10)*pH
@@ -697,7 +697,7 @@
             v_residue = rv[v_key]
             w_residue = rv[w_key]
 
-            for v_instance in v_residue.instances.values():
+            for v_instance in list(v_residue.instances.values()):
                 w_instances = list(w_residue.instances.values())
                 min_energy = min(self.normalized_interaction_energies[v_instance, w_instance]
                                  for w_instance in w_instances)
@@ -709,11 +709,11 @@
                         self.normalized_interaction_energies[w_instance, v_instance] -= min_energy
 
         #Normalize the instance energies.
-        for residue in self.residue_variables.values():
-            min_energy = min(instance.energyNF for instance in residue.instances.values())
+        for residue in list(self.residue_variables.values()):
+            min_energy = min(instance.energyNF for instance in list(residue.instances.values()))
             self.normalized_constant_energy += min_energy
             if min_energy != sys.float_info.max:
-                for instance in residue.instances.values():
+                for instance in list(residue.instances.values()):
                     instance.energyNF -= min_energy
 
 
@@ -722,11 +722,11 @@
         Used primarily for testing."""
         self.interaction_energies_for_ph = self.interaction_energies.copy()
 
-        for residue in self.residue_variables.values():
+        for residue in list(self.residue_variables.values()):
             residue.instances["DEPROTONATED"].energy_with_ph = residue.instances["DEPROTONATED"].energy
             residue.instances["PROTONATED"].energy_with_ph = residue.instances["PROTONATED"].energy + math.log(10)*pH
 
-        for v, w in combinations(iter(self.residue_variables.items()), 2):
+        for v, w in combinations(iter(list(self.residue_variables.items())), 2):
             v_key, v_residue = v
             w_key, w_residue = w
             is_same_his = (v_key[1:] == w_key[1:] and v_key[0] in ("HId", "HIe") and w_key[0] in ("HId", "HIe"))
--- a/pdb2pka/graph_cut/titration_curve.py
+++ b/pdb2pka/graph_cut/titration_curve.py
@@ -1,6 +1,6 @@
-from __future__ import print_function
-from graph import ProteinGraph
-from uncertainty import resolve_uncertainty
+
+from .graph import ProteinGraph
+from .uncertainty import resolve_uncertainty
 from collections import defaultdict
 import math
 from pprint import pprint
@@ -23,12 +23,12 @@
        normal_form - Dump the normal form part of the state."""
     rv = pc.residue_variables
     ie = pc.normalized_interaction_energies if normal_form else pc.interaction_energies_for_ph
-    for v_residue in rv.values():
-        for v_instance in v_residue.instances.values():
-            for w_residue in rv.values():
+    for v_residue in list(rv.values()):
+        for v_instance in list(v_residue.instances.values()):
+            for w_residue in list(rv.values()):
                 if v_residue == w_residue:
                     continue
-                for w_instance in w_residue.instances.values():
+                for w_instance in list(w_residue.instances.values()):
                     out_file.write(str((v_instance, w_instance)) + " " + str(round(ie[v_instance, w_instance],4)) + '\n')
     keys = list(pc.residue_variables.keys())
     for key in keys:
@@ -127,7 +127,7 @@
         new_labeling = resolve_uncertainty(protein_complex, labeling, uncertain, verbose=True)
 
         curve_values = get_curve_values(protein_complex, new_labeling, pH)
-        for key, value in curve_values.items():
+        for key, value in list(curve_values.items()):
             curves[key].append((pH, value))
 
     return curves
@@ -140,7 +140,7 @@
 
     aH = math.pow(10, -pH)
 
-    for key, residue in protein_complex.residue_variables.items():
+    for key, residue in list(protein_complex.residue_variables.items()):
         name, chain, location = key
 
         if name in ("HId", "HIe"):
--- a/pdb2pka/graph_cut/uncertainty.py
+++ b/pdb2pka/graph_cut/uncertainty.py
@@ -1,4 +1,4 @@
-from __future__ import print_function
+
 import sys
 from itertools import product
 import random
--- a/pdb2pka/graph_cut/utils.py
+++ b/pdb2pka/graph_cut/utils.py
@@ -18,8 +18,8 @@
         if flipped_inter_avg is not None:
             diff = abs(inter_avg - flipped_inter_avg)
             if diff > 0.0:
-                print group1_type, group1_chain, group1_loc, group1_state
-                print group2_type, group2_chain, group2_loc, group2_state
+                print(group1_type, group1_chain, group1_loc, group1_state)
+                print(group2_type, group2_chain, group2_loc, group2_state)
 
 
 
@@ -50,10 +50,10 @@
     res_type = pKa.pKaGroup.name
     chain = pKa.residue.chainID
     location = str(pKa.residue.resSeq)
-    for state, energy in pKa.desolvation.iteritems():
+    for state, energy in pKa.desolvation.items():
         _process_desolv_or_background_line(protein_complex, res_type, chain, location, state, energy)
 
-    for state, energy in pKa.background.iteritems():
+    for state, energy in pKa.background.items():
         _process_desolv_or_background_line(protein_complex, res_type, chain, location, state, energy)
 
 def _process_desolv_or_background_line(protein_complex, res_type, chain, location, state_name, energy):
--- a/pdb2pka/inputgen_pKa.py
+++ b/pdb2pka/inputgen_pKa.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 #
 # $Id$
 #
@@ -57,7 +57,7 @@
         #
         # Set all the attributes
         #
-        for prop in defaults.keys():
+        for prop in list(defaults.keys()):
             setattr(self,prop,defaults[prop])
         return
 
@@ -87,7 +87,7 @@
         #
         minmax=[[999.9,-999.9],[999.9,-999.9],[999.9,-999.9]]
         for coord in coords:
-            for axis in xrange(3):
+            for axis in range(3):
                 if coord[axis]<minmax[axis][0]:
                     minmax[axis][0]=coord[axis]
                 if coord[axis]>minmax[axis][1]:
@@ -97,7 +97,7 @@
         #
         center=[0,0,0]
         extent=[0,0,0]
-        for axis in xrange(3):
+        for axis in range(3):
             extent[axis]=minmax[axis][1]-minmax[axis][0]
             center[axis]=extent[axis]/2.0+minmax[axis][0]
         return center,extent
@@ -138,7 +138,7 @@
             #
             # Not a known type
             #
-            raise 'unknown type',type
+            raise Exception('unknown type')
         return
         
     def getText_sub(self):
--- a/pdb2pka/ligand_topology.py
+++ b/pdb2pka/ligand_topology.py
@@ -9,7 +9,7 @@
     import numpy as Numeric
     
 from sets import Set
-from ligandclean.trial_templates import *
+from .ligandclean.trial_templates import *
 from types import *
 
 def length(vector):
@@ -47,7 +47,7 @@
             # Get the likely types from names
             #
             trivial_types=['N','O','C','H']
-            for atom in self.atoms.keys():
+            for atom in list(self.atoms.keys()):
                 if atom[0] in trivial_types:
                     if atom[0]!='H': # Get rid of all the hydrogens
                         self.atoms[atom]['type']=atom[0]
@@ -57,9 +57,9 @@
             # First approximation: Anything closer than 2.0 A is bonded
             #
             self.dists={}
-            for atom1 in self.atoms.keys():
+            for atom1 in list(self.atoms.keys()):
                 self.dists[atom1]={}
-                for atom2 in self.atoms.keys():
+                for atom2 in list(self.atoms.keys()):
                     if atom1==atom2:
                         continue
                     #
@@ -75,7 +75,7 @@
             #
             # Get the torsion angles
             #
-            atoms=self.atoms.keys()
+            atoms=list(self.atoms.keys())
             atoms.sort()
             for atom in atoms:
                 self.atoms[atom]['torsions']=self.get_torsions(atom)
@@ -124,7 +124,7 @@
             #
             # Get the torsion angles
             #
-            atoms=self.atoms.keys()
+            atoms=list(self.atoms.keys())
             atoms.sort()
             for atom in atoms:
                 self.atoms[atom]['torsions']=self.get_torsions(atom)
@@ -143,7 +143,7 @@
         # + determine their likely order (e.g. single, double, or triple)
         #
         ambs={}
-        atoms=self.atoms.keys()
+        atoms=list(self.atoms.keys())
         for atom_name in atoms:
             bonds=self.atoms[atom_name]['bonds']
             atype=self.atoms[atom_name]['type']
@@ -166,13 +166,13 @@
         #
         valences={'C':4,'O':2,'N':3}
         for atom in self.atoms:
-            print atom, self.atoms[atom]
+            print(atom, self.atoms[atom])
         #
         # ok, now it gets hairy
         #
-        print
-        print 'Guessing sybyl atom types'
-        for atom in self.atoms.keys():
+        print()
+        print('Guessing sybyl atom types')
+        for atom in list(self.atoms.keys()):
             stype=None
             at=self.atoms[atom]
             #
@@ -212,7 +212,7 @@
         # Do some postchecks
         # - right now only for Carboxylic acids
         #
-        for atom in self.atoms.keys():
+        for atom in list(self.atoms.keys()):
             at=self.atoms[atom]
             if at['sybylType']=='C.2':
                 #
@@ -233,11 +233,11 @@
         #
         # All Done
         #
-        atoms=self.atoms.keys()
+        atoms=list(self.atoms.keys())
         atoms.sort()
-        print '\nFinal Sybyl type results'
+        print('\nFinal Sybyl type results')
         for atom in atoms:
-            print atom,self.atoms[atom]['sybylType']
+            print(atom,self.atoms[atom]['sybylType'])
         return
 
     #
@@ -267,7 +267,7 @@
         tps.sort() # To agree with dictionary layout
         best_fit=2.00
         best_type=None
-        for bond in bond_props.keys():
+        for bond in list(bond_props.keys()):
             if bond[0]==tps[0] and bond[-1]==tps[1]:
                 if abs(dist-bond_props[bond][0])<best_fit:
                     best_fit=abs(dist-bond_props[bond][0])
@@ -303,7 +303,7 @@
         # Filter the torsions
         #
         # 
-        print 'Jens has to write the stuff for filtering torsions..\n'
+        print('Jens has to write the stuff for filtering torsions..\n')
         return possible_torsions
                         
     #
@@ -315,7 +315,7 @@
         # Make the lines for the pdb2pqr definition
         #
         self.numbers={}
-        atoms=self.atoms.keys()
+        atoms=list(self.atoms.keys())
         atoms.sort()
         number=0
         for atom in atoms:
@@ -423,12 +423,12 @@
         #
         # Look for simple substructures that would be titratable groups in the ligand
         #
-        atoms=self.atoms.keys()
+        atoms=list(self.atoms.keys())
         #
         # ring detection (including deleting redundancies & sorting issues)
         ring_list = []
         tmp=[]
-        for atom in self.atoms.keys():
+        for atom in list(self.atoms.keys()):
             temp_ring_list = []
             tmp.append(self.ring_detection(atom))
         #
@@ -452,12 +452,12 @@
                     del sorted_ring_list[i]
                 else:
                     last = sorted_ring_list[i]
-        print "# overall rings (including potentially fused rings) :", len(sorted_ring_list)
+        print("# overall rings (including potentially fused rings) :", len(sorted_ring_list))
         stop ## PC 03.01.06
         #
         #
         # assigning ring attribute for every ring atom
-        for at in self.atoms.keys():
+        for at in list(self.atoms.keys()):
             self.atoms[at]['in_ring'] = 0
         for rring in sorted_ring_list:
             for current_atom in rring:
@@ -480,7 +480,7 @@
                       at['ring_list'] = [rring]
                 elif rring not in at['ring_list']:
                       at['ring_list'].append(rring)
-        print  "# non-fused rings                                   :", non_fused_counter
+        print("# non-fused rings                                   :", non_fused_counter)
 
  
         
@@ -490,7 +490,7 @@
         match_list=[]
             #
         # match ligand atom type with atom type from template
-        for at in t.keys():
+        for at in list(t.keys()):
             if t[at]['sybylType'] == self.atoms[atom2match]['sybylType']:
                 match_list.append(at)
             if len(match_list) != 0:
@@ -515,14 +515,14 @@
                             # Now match simultaneously atom_type and neighbouring atom_types for ligand AND template
                             if len(ligand_set.difference(template_set)) == 0 and len(ligand_list) == len(nbs_in_template):
                                 for entry in matched_atom_in_template:
-                                    print "%3d"%(counter),"  Ligand %4s %5s %28s " \
+                                    print("%3d"%(counter),"  Ligand %4s %5s %28s " \
                                           %(at_lig,self.atoms[at_lig]['sybylType'],ligand_list),\
                                           "template %s %s %s %s" \
-                                          %(matched_atom_in_template,t[entry]['sybylType'],nbs_in_template,t[entry]['neighbours'])
+                                          %(matched_atom_in_template,t[entry]['sybylType'],nbs_in_template,t[entry]['neighbours']))
                                     for neighboured_template_atoms in t[entry]['neighbours']:
-                                        print neighboured_template_atoms,t[neighboured_template_atoms]['sybylType'],t[neighboured_template_atoms]['sybyl_neighbours']
+                                        print(neighboured_template_atoms,t[neighboured_template_atoms]['sybylType'],t[neighboured_template_atoms]['sybyl_neighbours'])
                                     for neighboured_ligand_atoms in self.atoms[at_lig]['lBondedAtoms']:
-                                        print neighboured_ligand_atoms.name, neighboured_ligand_atoms.sybylType,neighboured_ligand_atoms.lBondedAtoms
+                                        print(neighboured_ligand_atoms.name, neighboured_ligand_atoms.sybylType,neighboured_ligand_atoms.lBondedAtoms)
                                     stop
                 counter += 1
 
@@ -530,15 +530,15 @@
             #
             # match ligand atom type with atom type from template
             if atom2match == "F14":
-                print "YYY_atom2match_YYY", atom2match
+                print("YYY_atom2match_YYY", atom2match)
 #            print "alrvis",len(already_visited),already_visited
             if matching_template == {}:
                 matching_template['MatchedFragments'] = {}
             if len(stored_nbs_of_atom2match) != 0 and stored_nbs_of_atom2match[-1] == atom2match:
-                print "bis zum erbrechen schreien!!!!", self.atoms[atom2match]['bonds']
+                print("bis zum erbrechen schreien!!!!", self.atoms[atom2match]['bonds'])
                 for e in  self.atoms[atom2match]['bonds']:
                     atom2match = e
-            for at in t.keys():
+            for at in list(t.keys()):
                 # TODO:matching ALL atom types in template => gives a match_list
                 if t[at]['sybylType'] == self.atoms[atom2match]['sybylType'] \
                    and atom2match not in already_visited:
@@ -587,9 +587,9 @@
                         else:
                             matched_atom_types2(nbat,t)
                 else:
-                    print "sybylType s don't match", atom2match
+                    print("sybylType s don't match", atom2match)
             # 2nd loop to go over to the neighboured atoms
-            for at in t.keys():
+            for at in list(t.keys()):
                 if atom2match in already_visited:
                     for nbat in self.atoms[atom2match]['bonds']:
                         if nbat in already_visited:
@@ -607,7 +607,7 @@
 
                         else:
                             matched_atom_types2(nbat,t)
-            print "\t\t\tlen alrvis %3d" % (len(already_visited))
+            print("\t\t\tlen alrvis %3d" % (len(already_visited)))
 
         def createsybyllistonthefly(lig_atom):
             # look in matched_atom_types2 - line 656
@@ -623,7 +623,7 @@
                 if ent_lig not in already_visited:
                     already_visited.append(ent_lig)
                     # look for matching neighbours
-                    for at in t.keys():
+                    for at in list(t.keys()):
                         if t[at]['sybylType'] == self.atoms[ent_lig]['sybylType'] and ent_lig not in matchlist:
                             matchlist.append(at)
                     for matches in matchlist:
@@ -634,7 +634,7 @@
                                 if putative_next_atom2match not in putative_next_a2m_list:
                                     putative_next_a2m_list.append(putative_next_atom2match)
                         else:
-                             print "sybyl neighbours don't match"
+                             print("sybyl neighbours don't match")
                 else:
                     # what's here?
                     pass
@@ -645,7 +645,7 @@
 
         def matchatomtypeintemplateandgetliglist(atom2match,t,stored_nbs_of_atom2match=[],been_here_flag=False,\
                                                  already_visited=[],hit_list=[]):
-            print atom2match,"hit_list",hit_list,been_here_flag
+            print(atom2match,"hit_list",hit_list,been_here_flag)
             putative_next_a2m_list = []
             # we don't want to miss the nbs of a matched atom (see [1])
             if atom2match in stored_nbs_of_atom2match:
@@ -653,15 +653,15 @@
                  gothroughallnbsofmatchlistatom(stored_nbs_of_atom2match,t,already_visited,hit_list)
             # does this really work? - to which position of the routine do we go now?
             if been_here_flag == True:
-                print "it's true...", putative_next_a2m_list
+                print("it's true...", putative_next_a2m_list)
                 for next_at in putative_next_a2m_list:
-                    print "TRUE (been_here_flag)", putative_next_a2m_list
+                    print("TRUE (been_here_flag)", putative_next_a2m_list)
                     matchatomtypeintemplateandgetliglist(next_at,t,been_here_flag=True)
             matchlist = []
-            for at in t.keys():
+            for at in list(t.keys()):
                 if t[at]['sybylType'] == self.atoms[atom2match]['sybylType'] and atom2match not in already_visited:
                     already_visited.append(atom2match)
-                    print "we found a match for %4s " %(atom2match)
+                    print("we found a match for %4s " %(atom2match))
                     matchlist.append(at)
             # look for sybylnbs of all stored entries in matchlist
             for entries in matchlist:
@@ -669,7 +669,7 @@
                 if len(Set(t[entries]['sybyl_neighbours']).difference(Set(createsybyllistonthefly(atom2match)))) == 0:
                     hit_list.append(atom2match)
                     stored_nbs_of_atom2match = self.atoms[atom2match]['bonds']
-                    print "nbs %s of hit %s" %(stored_nbs_of_atom2match,atom2match)
+                    print("nbs %s of hit %s" %(stored_nbs_of_atom2match,atom2match))
                     for nbs in stored_nbs_of_atom2match:
                         # call itself!
                         #
@@ -682,7 +682,7 @@
             matchatomtypeintemplateandgetliglist(start_atom,t)
 #            matched_atom_types2(start_atom,t)
 
-        for current_template in templates.keys():
+        for current_template in list(templates.keys()):
             match2(templates[current_template],atoms,start_atom=atoms[4])  # start_atom should be 0
 #            match(templates[current_template],atoms)
         
--- a/pdb2pka/ligandclean/ligff.py
+++ b/pdb2pka/ligandclean/ligff.py
@@ -1,5 +1,5 @@
 from src.forcefield import *
-from peoe_PDB2PQR import PEOE as calc_charges
+from .peoe_PDB2PQR import PEOE as calc_charges
 from src.pdb import *
 from src.definitions import *
 from pdb2pka import NEWligand_topology
@@ -86,7 +86,7 @@
     ligand_titratable_groups=X.find_titratable_groups()
 
     if verbose:
-        print "ligand_titratable_groups", ligand_titratable_groups
+        print("ligand_titratable_groups", ligand_titratable_groups)
     #
     # Append the ligand data to the end of the PDB data
     #
@@ -166,7 +166,7 @@
         elif ff == "parse":
             defpath = PARSE_FILE
         else:
-            raise ValueError, "Invalid forcefield %s!" % ff
+            raise ValueError("Invalid forcefield %s!" % ff)
 
         if not os.path.isfile(defpath):
             for path in sys.path:
@@ -175,7 +175,7 @@
                     defpath = testpath
                     break
         if not os.path.isfile(defpath):
-            raise ValueError, "Unable to find forcefield %s!" % defpath
+            raise ValueError("Unable to find forcefield %s!" % defpath)
 
         file = open(defpath, 'rU')
         lines = file.readlines()
@@ -268,7 +268,7 @@
                    "Cl": 1.75}
 
 #Add lower case keys to more lenient matching.
-ParseRadiiDictLower = dict((key.lower(), value) for key, value in ParseRadiiDict.items()) 
+ParseRadiiDictLower = dict((key.lower(), value) for key, value in list(ParseRadiiDict.items())) 
 ParseRadiiDict.update(ParseRadiiDictLower)
 
 class ligand_charge_handler(MOL2MOLECULE):
@@ -289,7 +289,7 @@
             #print "newly_calced - net charge %1.4f" %(qqqgesges)
             #print
         else:
-            atoms_last_calc=self.ligand_props.keys()
+            atoms_last_calc=list(self.ligand_props.keys())
             #
             # Get the atoms presently in the pdb2pqr array
             #
@@ -317,30 +317,30 @@
                 #
                 for atom in atoms_now:
                     if not atom in atoms_last_calc:
-                        print 'This atom was missing before',atom
-                        print 'If it is a hydrogen for the titratable, we need to create a bond entry!'
+                        print('This atom was missing before',atom)
+                        print('If it is a hydrogen for the titratable, we need to create a bond entry!')
                         # We should be here only if is a titratable
                         for current_atom in atoms_now:
                             # check if it't a titratable H
                             for res_atoms in residue.atoms:
                                 if current_atom == res_atoms.name and "titratableH"  in dir(res_atoms):
-                                    print "nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn"
-                                    print "been here"
+                                    print("nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn")
+                                    print("been here")
                                     for ResAtoms in residue.atoms:
                                         ResAtoms.formalcharge = 0.0
                                     self.recalc_charges(residue)
                 for atom in atoms_last_calc:
                     if not atom in atoms_now:
-                        print 'This atom used to be here, but is now missing',atom
+                        print('This atom used to be here, but is now missing',atom)
                 #self.recalc_charges(residue)
                 xxxnetq = 0.0
                 for xxx in residue.atoms:
-                    print "after neutralizing %s  %1.4f" %(xxx.name, xxx.charge)
+                    print("after neutralizing %s  %1.4f" %(xxx.name, xxx.charge))
                     xxxnetq = xxxnetq+xxx.charge
-                print '-----------------------'
-                print "net charge: %1.4f" % (xxxnetq)
-                print
-                print
+                print('-----------------------')
+                print("net charge: %1.4f" % (xxxnetq))
+                print()
+                print()
             else:
                 # Yes - nothing to do
                 pass
@@ -378,10 +378,10 @@
         for at in residue.atoms: # WAS: self.lAtoms:
             ele = at.sybylType.split('.')[0]
             charge = at.charge
-            if ParseRadiiDict.has_key(ele):
+            if ele in ParseRadiiDict:
                 radius = ParseRadiiDict[ele]
             else:
-                raise 'Please check ParseRadiiDict in ligff.py -- radius data not found for',ele
+                raise Exception('Please check ParseRadiiDict in ligff.py -- radius data not found for %s' % ele)
             #
             # Store the radii and charges
             #
--- a/pdb2pka/ligandclean/peoe_PDB2PQR.py
+++ b/pdb2pka/ligandclean/peoe_PDB2PQR.py
@@ -67,7 +67,7 @@
 }
 
 #Add lower case keys to more lenient matching.
-ChargetermsLower = dict((key.lower(), value) for key, value in Chargeterms.items()) 
+ChargetermsLower = dict((key.lower(), value) for key, value in list(Chargeterms.items())) 
 Chargeterms.update(ChargetermsLower)
 
 def PEOE( atoms, damp=0.778, k=1.56):
@@ -87,9 +87,9 @@
     abs_qges = 0.0
     counter = 0
     for a in atoms.atoms:
-    	sybylType = a.sybylType.lower()
-        if not Chargeterms.has_key(sybylType):
-            raise KeyError, 'PEOE Error: Atomtype <%s> not known, treating atom %s as dummy' % (a.sybylType, a.name)
+        sybylType = a.sybylType.lower()
+        if sybylType not in Chargeterms:
+            raise KeyError('PEOE Error: Atomtype <%s> not known, treating atom %s as dummy' % (a.sybylType, a.name))
         if a.sybylType == 'O.3':
             a.chi   = Chargeterms['O.OH'][0]
             a.abc   = Chargeterms['O.OH']
--- a/pdb2pka/pKaIO_compat.py
+++ b/pdb2pka/pKaIO_compat.py
@@ -12,7 +12,7 @@
 #try:
 import os
 
-from pKa_utility_functions_compat import *
+from .pKa_utility_functions_compat import *
 
 class pKaIO:
 
@@ -95,10 +95,10 @@
             filename=self.pkafile
         import string
         if not filename:
-            print filename,'is not a filename'
+            print(filename,'is not a filename')
             os._exit(0)
         if not os.path.isfile(filename):
-            raise 'File does not exist:',filename
+            raise Exception('File does not exist: %s' % filename)
         fd=open(filename)
         lines=fd.readlines()
         fd.close()
@@ -112,7 +112,7 @@
         else:
             raise ValueError('Unknown format')
         if string.lower(string.strip(lines[1]))!=string.lower('Format 1.0'):
-            raise 'unknown format: ',lines[1]
+            raise Exception('unknown format: %s' % lines[1])
         # Next line is text
         linenumber=3
         pKa={}
@@ -246,13 +246,13 @@
         fd.write('%s pKa File\n' %format)
         fd.write('Format 1.0\n')
         fd.write('      Group            pKa value  Model pK    dpK     dDesolv    dBack   dElec\n')
-        groups=data.keys()
+        groups=list(data.keys())
         # Sort accroding to the residue sequence number
         newgroup ={}
         for g in groups:
             newg = g.split('_')[2]
             newgroup[int(newg),g]=g
-        newgroupkeys = newgroup.keys()
+        newgroupkeys = list(newgroup.keys())
         newgroupkeys.sort()
         groups = []
         for k in newgroupkeys:
@@ -263,7 +263,7 @@
         # ---------
         #
         for group in groups:
-            if data.has_key(group):
+            if group in data:
                 this_data=data[group]
                 fd.write('%15s      %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  \n' %(self.WI_res_text(group,format),
                                                                                     this_data['pKa'],
@@ -295,7 +295,7 @@
             filename=self.titcurvfile
         import string
         if not os.path.isfile(filename):
-            raise 'File does not exist:',filename
+            raise Exception('File does not exist: %s' % filename)
         fd=open(filename)
         lines=fd.readlines()
         fd.close()
@@ -307,9 +307,9 @@
         elif string.lower(string.strip(lines[0]))==string.lower('pdb2pka Titration Curve File'):
             format='pdb2pka'
         else:
-            raise 'Unknown format'
+            raise Exception('Unknown format')
         if string.lower(string.strip(lines[1]))!=string.lower('Format 1.0'):
-            raise 'unknown format: ',lines[1]
+            raise Exception('unknown format: %s' % lines[1])
         phvals=string.split(lines[2])
         phstart=string.atof(phvals[0])
         phend=string.atof(phvals[1])
@@ -361,27 +361,27 @@
         #
         # Extract some data from the dictionary
         #
-        residues=data.keys()
-        phvals=data[residues[0]].keys()
+        residues=list(data.keys())
+        phvals=list(data[residues[0]].keys())
         phvals.sort()
         for residue in residues:
-            newpHvals=data[residue].keys()
+            newpHvals=list(data[residue].keys())
             newpHvals.sort()
             if newpHvals!=phvals:
-                print phvals
-                print newpHvals
-                raise 'Dictionary does not contain identical pH values'
+                print(phvals)
+                print(newpHvals)
+                raise Exception('Dictionary does not contain identical pH values')
         #
         # Check that a pKa value is in the pH values
         #
         for residue in residues:
-            if not data[residue].has_key('pKa'):
+            if 'pKa' not in data[residue]:
                 #print 'No pKa value found. Setting to zero!! - Jens change this!!'
                 data[residue]['pKa']=0.0
         #
         # Find the pH-start, stop and step
         #
-        phvals=data[residues[0]].keys()
+        phvals=list(data[residues[0]].keys())
         phvals.sort()
         phstart=phvals[0]
         phstop=phvals[-2]
@@ -397,13 +397,13 @@
         # Start pH, end pH, pH step
         #
         fd.write('%6.3f %7.3f %6.3f\n' %(phstart,phstop,phstep))
-        residues=data.keys()
+        residues=list(data.keys())
         # Sort accroding to the residue sequence number
         newresidue ={}
         for r in residues:
             newr = r.split('_')[2]
             newresidue[int(newr),r]=r
-        newresiduekeys = newresidue.keys()
+        newresiduekeys = list(newresidue.keys())
         newresiduekeys.sort()
         residues = []
         for k in newresiduekeys:
@@ -439,11 +439,11 @@
             if self.matrix_file:
                 filename=self.matrix_file
             else:
-                raise 'No matrix filename given'
+                raise Exception('No matrix filename given')
         #
         import os, string
         if not os.path.isfile(filename):
-            raise "File not found",filename
+            raise Exception("File not found %s" % filename)
         fd=open(filename)
         lines=fd.readlines()
         fd.close()
@@ -459,9 +459,9 @@
         elif string.lower(string.strip(lines[0]))==string.lower('pdb2pka Interaction Matrix File'):
             format='WHAT IF'
         else:
-            raise 'Unknown format'
+            raise Exception('Unknown format')
         if not string.strip(lines[1])=='Format 1.0':
-            raise 'Wrong format',lines[1]
+            raise Exception('Wrong format %s' % lines[1])
         x=1
         done=None
         partners=None
@@ -500,7 +500,7 @@
                     partners=np
                 else:
                     if partners!=np:
-                        raise 'Number of partners changes:',np
+                        raise Exception('Number of partners changes: %s' % np)
                 self.matrix[resid]={}
                 #
                 # Now read all the interactions with the partners
@@ -577,7 +577,7 @@
         fd=open(filename,'w')
         fd.write('%s Interaction Matrix File\n' %format)
         fd.write('Format 1.0\n')
-        groups=self.matrix.keys()
+        groups=list(self.matrix.keys())
         groups.sort()
 
         num_groups=len(groups)
@@ -594,7 +594,7 @@
         for g in newgroups:
             newg = g.split('_')[2]
             newnewgroup[int(newg),g]=g
-        newnewgroupkeys = newnewgroup.keys()
+        newnewgroupkeys = list(newnewgroup.keys())
         newnewgroupkeys.sort()
         newgroups = []
         for k in newnewgroupkeys:
@@ -605,19 +605,19 @@
         # ---------
         #
         for group in newgroups:
-            if self.matrix.has_key(group):
+            if group in self.matrix:
                 fd.write('%s      %7.4f\n' %(self.WI_res_text(group,format),float(num_groups)))
                 self.write_section(group,fd,format)
                 written[group]=1
                 #
                 # Is there a terminal group associated with this residue?
                 #
-                if self.matrix.has_key('T'+group):
+                if 'T'+group in self.matrix:
                     fd.write('%s      %7.4f\n' %(self.WI_res_text('T'+group),float(num_groups)))
                     self.write_section('T'+group,fd,format)
                     written['T'+group]=1
             else:
-                if self.matrix.has_key('T'+group) and not written.has_key('T'+group):
+                if 'T'+group in self.matrix and 'T'+group not in written:
                     fd.write('%s      %7.4f\n' %(self.WI_res_text('T'+group,format),float(num_groups)))
                     self.write_section('T'+group,fd,format)
                     written['T'+group]=1
@@ -632,19 +632,19 @@
         """Write an interaction energy matrix in pdb2pka format"""
         """At the moment, we just reformat and write a WHAT IF file"""
         self.matrix={}
-        for group1 in matrix.keys():
+        for group1 in list(matrix.keys()):
             self.matrix[group1.uniqueid]={}
-            for tit1 in matrix[group1].keys():
-                for state1 in matrix[group1][tit1].keys():
+            for tit1 in list(matrix[group1].keys()):
+                for state1 in list(matrix[group1][tit1].keys()):
                     sub_m=matrix[group1][tit1][state1]
-                    for group2 in sub_m.keys():
-                        if not self.matrix[group1.uniqueid].has_key(group2.uniqueid):
+                    for group2 in list(sub_m.keys()):
+                        if group2.uniqueid not in self.matrix[group1.uniqueid]:
                             self.matrix[group1.uniqueid][group2.uniqueid]=[]
-                        for tit2 in sub_m[group2].keys():
-                            for state2 in sub_m[group2][tit2].keys():
+                        for tit2 in list(sub_m[group2].keys()):
+                            for state2 in list(sub_m[group2][tit2].keys()):
                                 self.matrix[group1.uniqueid][group2.uniqueid].append(sub_m[group2][tit2][state2])
-        for group1 in self.matrix.keys():
-            for group2 in self.matrix[group1].keys():
+        for group1 in list(self.matrix.keys()):
+            for group2 in list(self.matrix[group1].keys()):
                 sum=0.0
                 for val in self.matrix[group1][group2]:
                     sum=sum+val
@@ -657,7 +657,7 @@
     #
 
     def write_section(self,group,fd,format):
-        groups_tmp=self.matrix[group].keys()
+        groups_tmp=list(self.matrix[group].keys())
         groups2=[]
         for group_x in groups_tmp:
             if group_x[0]=='T':
@@ -670,7 +670,7 @@
         for g in groups2:
             newg = g.split('_')[2]
             newgroup[int(newg),g]=g
-        newgroupkeys = newgroup.keys()
+        newgroupkeys = list(newgroup.keys())
         newgroupkeys.sort()
         groups2 = []
         for k in newgroupkeys:
@@ -678,19 +678,19 @@
 
         written={}
         for group2 in groups2:
-            if self.matrix[group].has_key(group2):
+            if group2 in self.matrix[group]:
                 fd.write('%s      %7.4f\n' %(self.WI_res_text(group2,format),self.matrix[group][group2][0]))
                 fd.write('%7.4f\n%7.4f\n%7.4f\n'%(self.matrix[group][group2][1],self.matrix[group][group2][2],self.matrix[group][group2][3]))
                 written[group2]=1
                 #
                 # Is there a terminal group associated with this residue?
                 #
-                if self.matrix[group].has_key('T'+group2):
+                if 'T'+group2 in self.matrix[group]:
                     fd.write('%s      %7.4f\n' %(self.WI_res_text('T'+group2,format),self.matrix[group]['T'+group2][0]))
                     fd.write('%7.4f\n%7.4f\n%7.4f\n'%(self.matrix[group]['T'+group2][1],self.matrix[group]['T'+group2][2],self.matrix[group]['T'+group2][3]))
                     written['T'+group2]=1
             else:
-                if self.matrix[group].has_key('T'+group2) and not written.has_key('T'+group2):
+                if 'T'+group2 in self.matrix[group] and 'T'+group2 not in written:
                     fd.write('%s      %7.4f\n' %(self.WI_res_text('T'+group2,format),self.matrix[group]['T'+group2][0]))
                     fd.write('%7.4f\n%7.4f\n%7.4f\n'%(self.matrix[group]['T'+group2][1],self.matrix[group]['T'+group2][2],self.matrix[group]['T'+group2][3]))
                     written['T'+group2]=1 
@@ -706,14 +706,14 @@
             if self.desolv_file:
                 filename=self.desolv_file
             else:
-                raise 'No desolv filename given'
+                raise Exception('No desolv filename given')
         #
         #
         # This subroutine reads a DESOLV file
         #
         import os, string
         if not os.path.isfile(filename):
-            raise "File not found",filename
+            raise Exception("File not found %s" % filename)
         fd=open(filename)
         lines=fd.readlines()
         fd.close()
@@ -729,7 +729,7 @@
         elif string.strip(lines[0])=='pdb2pka Desolvation Energy File':
             format='pdb2pka'
         else:
-            raise Exception,'Unknown format:'+string.strip(lines[0])
+            raise Exception('Unknown format:'+string.strip(lines[0]))
         #
         # Call the generic read routine
         #
@@ -748,11 +748,11 @@
             if self.backgr_file:
                 filename=self.backgr_file
             else:
-                raise 'No matrix filename given'
+                raise Exception('No matrix filename given')
         #
         import os, string
         if not os.path.isfile(filename):
-            raise "File not found",filename
+            raise Exception("File not found %s" % filename)
         fd=open(filename)
         lines=fd.readlines()
         fd.close()
@@ -768,7 +768,7 @@
         elif string.strip(lines[0])=='pdb2pka Background Energy File':
             format='pdb2pka'
         else:
-            raise Exception,'Unknown format:'+string.strip(lines[0])
+            raise Exception('Unknown format:'+string.strip(lines[0]))
         #
         # Call the generic read routine
         #
@@ -830,13 +830,13 @@
         fd=open(filename,'w')
         fd.write('%s Desolvation Energy File\n' %format)
         fd.write('Format 1.0\n')
-        groups=self.desolv.keys()
+        groups=list(self.desolv.keys())
         # Sort accroding to the residue sequence number
         newgroup ={}
         for g in groups:
             newg = g.split('_')[2]
             newgroup[int(newg),g]=g
-        newgroupkeys = newgroup.keys()
+        newgroupkeys = list(newgroup.keys())
         newgroupkeys.sort()
         groups = []
         for k in newgroupkeys:
@@ -846,7 +846,7 @@
         # ---------
         #
         for group in groups:
-            if self.desolv.has_key(group):
+            if group in self.desolv:
                 fd.write('%s      %7.4f\n' %(self.WI_res_text(group,format),float(self.desolv[group])))
                 written[group]=1
         fd.write('End of file\n')
@@ -863,13 +863,13 @@
         fd=open(filename,'w')
         fd.write('%s Background Energy File\n' %format)
         fd.write('Format 1.0\n')
-        groups=self.backgr.keys()
+        groups=list(self.backgr.keys())
         # Sort accroding to the residue sequence number
         newgroup ={}
         for g in groups:
             newg = g.split('_')[2]
             newgroup[int(newg),g]=g
-        newgroupkeys = newgroup.keys()
+        newgroupkeys = list(newgroup.keys())
         newgroupkeys.sort()
         groups = []
         for k in newgroupkeys:
@@ -879,7 +879,7 @@
         # ---------
         #
         for group in groups:
-            if self.backgr.has_key(group):
+            if group in self.backgr:
                 fd.write('%s      %7.4f\n' %(self.WI_res_text(group,format),float(self.backgr[group])))
                 written[group]=1
         fd.write('End of file\n')
@@ -955,4 +955,4 @@
                 line='TERMINAL GROUP:\n'+line
             return line
         else:
-            raise Exception,'Unknown format:'+format
+            raise Exception('Unknown format:'+format)
--- a/pdb2pka/pKa_utility_functions_compat.py
+++ b/pdb2pka/pKa_utility_functions_compat.py
@@ -163,7 +163,7 @@
 def istitratable(uniqueid):
     import string
     type=string.split(uniqueid,':')[-1]
-    if acidbase.has_key(type):
+    if type in acidbase:
         return 1
     else:
         return None
--- a/pdb2pka/pMC_mult.py
+++ b/pdb2pka/pMC_mult.py
@@ -25,7 +25,7 @@
     if (name == "thisown"): return self.this.own()
     method = class_type.__swig_getmethods__.get(name,None)
     if method: return method(self)
-    raise AttributeError,name
+    raise AttributeError(name)
 
 def _swig_repr(self):
     try: strthis = "proxy of " + self.this.__repr__()
@@ -34,7 +34,7 @@
 
 import types
 try:
-    _object = types.ObjectType
+    _object = object
     _newclass = 1
 except AttributeError:
     class _object : pass
@@ -47,7 +47,7 @@
     __setattr__ = lambda self, name, value: _swig_setattr(self, PySwigIterator, name, value)
     __swig_getmethods__ = {}
     __getattr__ = lambda self, name: _swig_getattr(self, PySwigIterator, name)
-    def __init__(self): raise AttributeError, "No constructor defined"
+    def __init__(self): raise AttributeError("No constructor defined")
     __repr__ = _swig_repr
     __swig_destroy__ = _pMC_mult.delete_PySwigIterator
     __del__ = lambda self : None;
--- a/pdb2pka/pka_help.py
+++ b/pdb2pka/pka_help.py
@@ -28,14 +28,14 @@
     #
     # Call our little C++ module
     #
-    import pMC_mult
+    from . import pMC_mult
     FAST=pMC_mult.MC(intpkas,linear,acidbase,state_counter,is_charged)
     FAST.set_MCsteps(int(mcsteps))
-    print 'Calculating intrinsic pKa value'
+    print('Calculating intrinsic pKa value')
     pKavals=FAST.calc_pKas(phstart,phend,phstep)
     count=0
     intpka=pKavals[0]
-    print 'Simulated intrinsic pKa value: %5.2f' %intpka
+    print('Simulated intrinsic pKa value: %5.2f' %intpka)
     count=1
     #
     # Get the charges
@@ -57,8 +57,8 @@
     if pKavals[count+1]==999.0 and pKavals[count+2]==-999.0:
         count=count+2
     else:
-        print 'Something is wrong'
-        print pKavals[count:count+30]
+        print('Something is wrong')
+        print(pKavals[count:count+30])
         raise Exception('Incorrect data format from pMC_mult')
     return intpka
 
@@ -73,7 +73,7 @@
                 ##check if the record is not a water in which case we will print a warning
                 from src.aa import WAT
                 if not record.resName in WAT.water_residue_names:
-                    print "Warning!: HETATM record %s " % record.resName + "%s that is not a water is being dropped\n  " % record.element
+                    print("Warning!: HETATM record %s " % record.resName + "%s that is not a water is being dropped\n  " % record.element)
                     ##raw_input("Press enter to continue...")
                 continue
             if isinstance(record, (ATOM, ANISOU, SIGUIJ, SIGATM)):
--- a/pdb2pka/pka_old.py
+++ b/pdb2pka/pka_old.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
 #
 # pKa calculations with APBS
 #
@@ -12,9 +12,9 @@
 debug=False
 import optparse
 import sys, os
-from pKa_base import *
+from .pKa_base import *
 
-print __file__
+print(__file__)
 import os
 try:
     file_name=__file__
@@ -55,7 +55,7 @@
 from src.routines import *
 from src.protein import *
 from src.server import *
-from StringIO import *
+from io import *
 from src.hydrogens import *
 
 
@@ -63,11 +63,11 @@
     #
     # Print usage guidelines
     #
-    print 'Usage: pka.py --ff <forcefield> --lig <ligand in MOL2> --pdie <protein diel cons> --maps <1 for using provided 3D maps; 2 for genereting new maps>'
-    print '--xdiel <xdiel maps> --ydiel <ydiel maps> --zdiel <zdiel maps> --kappa <ion-accessibility map> '
-    print '--smooth <st.dev [A] of Gaussian smooting of 3D maps at the boundary, bandthwith=3 st.dev> <pdbfile>'
-    print 'Force field can be amber, charmm and parse'
-    print
+    print('Usage: pka.py --ff <forcefield> --lig <ligand in MOL2> --pdie <protein diel cons> --maps <1 for using provided 3D maps; 2 for genereting new maps>')
+    print('--xdiel <xdiel maps> --ydiel <ydiel maps> --zdiel <zdiel maps> --kappa <ion-accessibility map> ')
+    print('--smooth <st.dev [A] of Gaussian smooting of 3D maps at the boundary, bandthwith=3 st.dev> <pdbfile>')
+    print('Force field can be amber, charmm and parse')
+    print()
     return
 
 #
@@ -83,9 +83,9 @@
             routines:   The routines object as generated by PDB2PQR
             forcefield: The forcefield object as generated by PDB2PQR
     """
-    print
-    print 'PDB2PQR pKa calculations'
-    print
+    print()
+    print('PDB2PQR pKa calculations')
+    print()
     
     import optparse
     parser = optparse.OptionParser()
@@ -252,7 +252,7 @@
     #
     # remove hydrogen atoms
     #
-    import pka_help
+    from . import pka_help
     pka_help.remove_hydrogens(pdbfilename)
     #
     # Get the PDBfile
@@ -268,11 +268,11 @@
 #         sys.exit(2)
 
     if len(errlist) != 0 and verbose:
-        print "Warning: %s is a non-standard PDB file.\n" %pdbfilename
-        print errlist
+        print("Warning: %s is a non-standard PDB file.\n" %pdbfilename)
+        print(errlist)
 
     if verbose:
-        print "Beginning PDB2PQR...\n"
+        print("Beginning PDB2PQR...\n")
     #
     # Read the definition file
     #
@@ -301,12 +301,12 @@
                 if os.path.isfile(ligand):
                     ligfd=open(ligand, 'rU')
                 else:
-                    print 'skipping ligand',ligand
+                    print('skipping ligand',ligand)
                     continue
                 #
                 # Read the ligand into Paul's code
                 #
-                from ligandclean import ligff
+                from .ligandclean import ligff
                 myProtein, myDefinition, Lig = ligff.initialize(myDefinition, ligfd, pdblist, verbose)
                 #
                 # Create the ligand definition from the mol2 data
@@ -387,9 +387,9 @@
     # Print something for some reason?
     #
     if verbose:
-        print "Created protein object -"
-        print "\tNumber of residues in protein: %s" % myProtein.numResidues()
-        print "\tNumber of atoms in protein   : %s" % myProtein.numAtoms()
+        print("Created protein object -")
+        print("\tNumber of residues in protein: %s" % myProtein.numResidues())
+        print("\tNumber of atoms in protein   : %s" % myProtein.numAtoms())
     #
     # Set up all other routines
     #
@@ -440,8 +440,8 @@
                 templist = []
                 Lig.make_up2date(residue)
                 net_charge=0.0
-                print 'Ligand',residue
-                print 'Atom\tCharge\tRadius'
+                print('Ligand',residue)
+                print('Atom\tCharge\tRadius')
                 for atom in residue.getAtoms():
                     if atom.mol2charge:
                         atom.ffcharge=atom.mol2charge
@@ -455,7 +455,7 @@
                     # Assign radius
                     #
                     atom.radius = Lig.ligand_props[atom.name]["radius"]
-                    print '%s\t%6.4f\t%6.4f' %(atom.name,atom.ffcharge,atom.radius)
+                    print('%s\t%6.4f\t%6.4f' %(atom.name,atom.ffcharge,atom.radius))
                     if atom in misslist:
                         misslist.pop(misslist.index(atom))
                         templist.append(atom)
@@ -480,7 +480,7 @@
                 #
                 # Print the net charge
                 #
-                print 'Net charge for ligand %s is: %5.3f' %(residue.name,net_charge)
+                print('Net charge for ligand %s is: %5.3f' %(residue.name,net_charge))
         #
         # Temporary fix; if ligand was successful, pull all ligands from misslist
        # Not sure if this is needed at all here ...? (Jens wrote this)
@@ -492,9 +492,9 @@
                 misslist.remove(atom)                    
     
     if verbose:
-        print "Created protein object (after processing myRoutines) -"
-        print "\tNumber of residues in protein: %s" % myProtein.numResidues()
-        print "\tNumber of atoms in protein   : %s" % myProtein.numAtoms()
+        print("Created protein object (after processing myRoutines) -")
+        print("\tNumber of residues in protein: %s" % myProtein.numResidues())
+        print("\tNumber of atoms in protein   : %s" % myProtein.numAtoms())
     #
     # Create the APBS input file
     #
@@ -503,17 +503,17 @@
 
     method=""
     split=0
-    import inputgen_pKa
+    from . import inputgen_pKa
     igen = inputgen_pKa.inputGen(pdbfilename)
     #
     # For convenience
     #
     igen.pdie = pdie
-    print 'Setting protein dielectric constant to ',igen.pdie
+    print( 'Setting protein dielectric constant to ',igen.pdie)
     igen.sdie=options.sdie
     igen.maps=maps
     if maps==1:
-        print "Using dielectric and mobile ion-accessibility function maps in PBE"
+        print( "Using dielectric and mobile ion-accessibility function maps in PBE")
         if xdiel:
             igen.xdiel = xdiel
         else:
@@ -533,7 +533,7 @@
             usage(2)
             sys.exit(0)
         
-        print 'Setting dielectric function maps: %s, %s, %s'%(igen.xdiel,igen.ydiel,igen.zdiel)
+        print('Setting dielectric function maps: %s, %s, %s'%(igen.xdiel,igen.ydiel,igen.zdiel))
         
         if kappa:
             igen.kappa = kappa
@@ -542,7 +542,7 @@
             usage(2)
             sys.exit(0)
             
-        print 'Setting mobile ion-accessibility function map to: ',igen.kappa
+        print('Setting mobile ion-accessibility function map to: ',igen.kappa)
         
         if sd:
             xdiel_smooth, ydiel_smooth, zdiel_smooth = smooth(xdiel,ydiel,zdiel)
@@ -561,7 +561,7 @@
 if __name__ == "__main__":
     
     (protein, routines, forcefield,apbs_setup, ligand_titratable_groups, maps, sd), options = startpKa()
-    import pka_routines
+    from . import pka_routines
     mypkaRoutines = pka_routines.pKaRoutines(protein, routines, forcefield, apbs_setup, maps, sd,
                                              pdbfile_name,
                                              options=options)
@@ -579,7 +579,7 @@
     #
     # What should we do?
     #
-    print 'Doing full pKa calculation'
+    print('Doing full pKa calculation')
     mypkaRoutines.runpKa()
 
 #     state=1
--- a/pdb2pka/pka_routines.py
+++ b/pdb2pka/pka_routines.py
@@ -10,21 +10,21 @@
 debug=False
 import os
 import sys
-import pKaIO_compat
-from pKa_base import *
-import cPickle
-import pMC_mult
+from . import pKaIO_compat
+from .pKa_base import *
+import pickle
+from . import pMC_mult
 import math
 import copy
 import string
 
-from graph_cut.utils import create_protein_complex_from_matrix, process_desolv_and_background, curve_for_one_group
-from graph_cut.titration_curve import get_titration_curves
-from graph_cut.create_titration_output import create_output
+from .graph_cut.utils import create_protein_complex_from_matrix, process_desolv_and_background, curve_for_one_group
+from .graph_cut.titration_curve import get_titration_curves
+from .graph_cut.create_titration_output import create_output
 
 if debug:
-    from Tkinter import *
-    from charge_mon import *
+    from tkinter import *
+    from .charge_mon import *
 
     CM=charge_mon()
 else:
@@ -32,7 +32,7 @@
 
 import shutil
 
-from pka_help import is_sameatom, titrate_one_group
+from .pka_help import is_sameatom, titrate_one_group
 from src.errors import PDB2PKAError
 
 #
@@ -56,7 +56,7 @@
 from src.hydrogens import hydrogenRoutines, hydrogenAmbiguity
 
 #import ligandclean.ligff
-from apbs import runAPBS
+from .apbs import runAPBS
 
 
 #
@@ -124,7 +124,7 @@
                 if os.path.isfile(self.state_files):
                     raise ValueError('Target directory is a file! Aborting.')
 
-                for output_file in self.output_files.values():
+                for output_file in list(self.output_files.values()):
                     if os.path.isfile(output_file):
                         os.remove(output_file)
 
@@ -161,7 +161,7 @@
     def insert_new_titratable_group(self,ligand_titratable_groups):
         """Insert the new titratable groups in to self.pkagroups"""
         group_type=ligand_titratable_groups['type']
-        if self.pKagroups.has_key(group_type):
+        if group_type in self.pKagroups:
             #
             # Now modify the group so that it will correspond to the group
             # we have in the ligand
@@ -190,7 +190,7 @@
                         #
                         # Change the name of the atom that the H is bound to
                         #
-                        if atom_map.has_key(conformation.boundatom):
+                        if conformation.boundatom in atom_map:
                             conformation.boundatom=atom_map[conformation.boundatom]
                         #
                         # Change the name of the hydrogen
@@ -201,7 +201,7 @@
                         # And then for the individual atom names
                         #
                         for atom in conformation.atoms:
-                            if atom_map.has_key(atom.name):
+                            if atom.name in atom_map:
                                 atom.name=atom_map[atom.name]
                             elif atom.name==oldhname:
                                 atom.name=conformation.hname
@@ -316,7 +316,7 @@
             ambiguity = pKa.amb
 
             titgroup='%s:%s:%s' %(residue.chainID, string.zfill(residue.resSeq,4),pKaGroup.name)
-            if not potentialDifference.has_key(titgroup):
+            if titgroup not in potentialDifference:
                 potentialDifference[titgroup]={}
                 for atom in self.protein.getAtoms():
                     if atom.name in ['N','H','C']:
@@ -402,12 +402,12 @@
         # Loop over each titration
         #
         self.all_potentials={}
-        if not self.matrix.has_key(pKa):
+        if pKa not in self.matrix:
             self.matrix[pKa]={}
             self.all_potentials[pKa]={}
         #
         for titration in pKaGroup.DefTitrations:
-            if not self.matrix[pKa].has_key(titration):
+            if titration not in self.matrix[pKa]:
                 self.matrix[pKa][titration]={}
                 self.all_potentials[pKa][titration]={}
             #
@@ -472,7 +472,7 @@
         #
         if os.path.isfile(intene_file_name):
             with open(intene_file_name) as fd:
-                savedict=cPickle.load(fd)
+                savedict=pickle.load(fd)
             #
             # pKD?
             #
@@ -480,7 +480,7 @@
                 try:
                     sys.stdout.flush()
                     with open(allpots_file_name) as fd:
-                        allsavedict=cPickle.load(fd)
+                        allsavedict=pickle.load(fd)
                     read_allpots=1
                 except EOFError:
                     self.routines.write('\n')
@@ -518,12 +518,12 @@
             #
             # Loop over each titration
             #
-            if not energies.has_key(pKa):
+            if pKa not in energies:
                 energies[pKa]={}
                 all_potentials[pKa]={}
             #
             for titration in pKaGroup.DefTitrations:
-                if not energies[pKa].has_key(titration):
+                if titration not in energies[pKa]:
                     energies[pKa][titration]={}
                     all_potentials[pKa][titration]={}
                 #
@@ -549,10 +549,10 @@
                     #
                     # Check if we have values for this calculation already
                     #
-                    if savedict.has_key(name):
+                    if name in savedict:
                         energies[pKa][titration][state]= savedict[name]
                         if mode=='pKD':
-                            if allsavedict.has_key(name):
+                            if name in allsavedict:
                                 all_potentials[pKa][titration][state]=allsavedict[name]
                         continue
                     #
@@ -675,10 +675,10 @@
         #
         if calculated_energy:
             with open(intene_file_name,'w') as fd:
-                cPickle.dump(savedict,fd)
+                pickle.dump(savedict,fd)
             if mode=='pKD' and not read_allpots:
                 with open(allpots_file_name,'w') as fd:
-                    cPickle.dump(allsavedict,fd)
+                    pickle.dump(allsavedict,fd)
         return energies,all_potentials
 
     #
@@ -728,8 +728,8 @@
         create_output(self.titcurves_dir, curves)
 
         pka_values, pH_values = self.find_pka_and_pH(curves)
-        print pka_values
-        print pH_values
+        print(pka_values)
+        print(pH_values)
         self.ph_at_0_5 = pH_values
 
         ln10=math.log(10)
@@ -812,7 +812,7 @@
         #
         X.desolv={}
         X.backgr={}
-        for name in pkas.keys():
+        for name in list(pkas.keys()):
             X.desolv[name]=pkas[name]['desolv']
             X.backgr[name]=pkas[name]['backgr']
 
@@ -842,7 +842,7 @@
 
         adjacent_data_points = 5
 
-        for name, curve in curves.iteritems():
+        for name, curve in curves.items():
             bad_curve = False
             curve_calc_point = 0.5
 
@@ -854,7 +854,7 @@
             #Check to see if we never cross 0.5 or -0.5 at all
             if all(charge_side) or not any(charge_side):
                 warning = "WARNING: UNABLE TO CACLCULATE PKA FOR {name}\n".format(name=name)
-                print warning,
+                print(warning, end=' ')
                 self.warnings.append(warning)
                 pKa_results[name] = pKa_value
                 continue
@@ -865,11 +865,11 @@
 
             if (True,False) not in side_pairs:
                 warning =  "WARNING: {name} DOES NOT EXHIBIT Henderson-Hasselbalch BEHAVIOR\n".format(name=name)
-                print warning,
+                print(warning, end=' ')
                 self.warnings.append(warning)
 
                 warning = "WARNING: {name} TITRATION CURVE IS BACKWARDS\n".format(name=name, calc=curve_calc_point)
-                print warning,
+                print(warning, end=' ')
                 self.warnings.append(warning)
                 cross_index = charge_side.index(True)
                 bad_curve = True
@@ -880,10 +880,10 @@
             #(False, True) means we've crossed back over the line and therefore our PKA value is in question.
             if (True,False) in side_pairs and (False,True) in side_pairs:
                 warning =  "WARNING: {name} DOES NOT EXHIBIT Henderson-Hasselbalch BEHAVIOR\n".format(name=name)
-                print warning,
+                print(warning, end=' ')
                 self.warnings.append(warning)
                 warning = "WARNING: {name} TITRATION CURVE CROSSES {calc} AT LEAST TWICE\n".format(name=name, calc=curve_calc_point)
-                print warning,
+                print(warning, end=' ')
                 self.warnings.append(warning)
 
                 bad_curve = True
@@ -899,11 +899,11 @@
                 pH_results[name] = ph_at_0_5
             except ZeroDivisionError:
                 warning = "WARNING: UNABLE TO CACLCULATE pH FOR {name}, Divide by zero.\n".format(name=name)
-                print warning,
+                print(warning, end=' ')
                 self.warnings.append(warning)
 
             if not bad_curve:
-                print "{name} exhibits Henderson-Hasselbalch behavior.".format(name=name)
+                print("{name} exhibits Henderson-Hasselbalch behavior.".format(name=name))
 
             #Calc pKa value
             start = max(0, prevous_cross_index-adjacent_data_points)
@@ -915,7 +915,7 @@
                 pKa_value = sum(pkas)/float(len(pkas))
             except ZeroDivisionError:
                 warning = "WARNING: UNABLE TO CACLCULATE PKA FOR {name}, Divide by zero.\n".format(name=name)
-                print warning,
+                print(warning, end=' ')
                 self.warnings.append(warning)
 
 
@@ -1191,7 +1191,7 @@
 #             intpka=titrate_one_group(name='%s' %(pKa.residue),intpkas=intpKas,is_charged=is_charged,acidbase=acidbase)
             curve = curve_for_one_group(pKa)
             pka_values, _ = self.find_pka_and_pH(curve)
-            intpka = pka_values.values()[0]
+            intpka = list(pka_values.values())[0]
             pKa.simulated_intrinsic_pKa=intpka
         return
 
@@ -1241,7 +1241,7 @@
         backgroundname = os.path.join(self.state_files,'background_interaction_energies.pickle')
         if os.path.isfile(backgroundname):
             with open(backgroundname) as fd:
-                savedict=cPickle.load(fd)
+                savedict=pickle.load(fd)
         else:
             savedict={}
 
@@ -1270,7 +1270,7 @@
                     # Set the name for this energy
                     #
                     name='%s_%s_%s_%s' %(titration.name,pKa.residue.chainID,pKa.residue.resSeq,self.get_state_name(titration.name,state))
-                    if savedict.has_key(name):
+                    if name in savedict:
                         pKa.background[self.get_state_name(titration.name,state)] = savedict[name]
                         continue
                     #
@@ -1416,10 +1416,10 @@
                     # Dump the pickle file
                     #
                     with open(backgroundname,'w') as fd:
-                        cPickle.dump(savedict,fd)
+                        pickle.dump(savedict,fd)
 
         with open(self.output_files['background_interaction_energies_file_path'] , 'w') as f:
-            keys = savedict.keys()
+            keys = list(savedict.keys())
             keys.sort()
             for key in keys:
                 value = savedict[key]
@@ -1442,7 +1442,7 @@
         desolvname=os.path.join(self.state_files, 'desolvation_energies.pickle')
         if os.path.isfile(desolvname):
             with open(desolvname) as fd:
-                savedict=cPickle.load(fd)
+                savedict=pickle.load(fd)
         else:
             savedict={}
         #
@@ -1503,7 +1503,7 @@
                     residue.stateboolean[self.get_state_name(titration.name, state)] = True
                     name='%s_%s_%s_%s' %(titration.name,pKa.residue.chainID,pKa.residue.resSeq,self.get_state_name(titration.name,state))
                     #
-                    if savedict.has_key(name):
+                    if name in savedict:
                         pKa.desolvation[self.get_state_name(titration.name,state)] = savedict[name]
                         continue
                     self.routines.write("---------> Calculating desolvation energy for residue %s state %s in solvent\n" %(residue.name,self.get_state_name(titration.name,state)))
@@ -1570,10 +1570,10 @@
                     # Dump a pickle file
                     #
                     with open(desolvname,'w') as fd:
-                        cPickle.dump(savedict,fd)
+                        pickle.dump(savedict,fd)
 
         with open(self.output_files['desolvation_energies_file_path'], 'w') as f:
-            keys = savedict.keys()
+            keys = list(savedict.keys())
             keys.sort()
             for key in keys:
                 value = savedict[key]
@@ -1803,7 +1803,7 @@
         #
         center={}
         extent={}
-        for axis in minmax.keys():
+        for axis in list(minmax.keys()):
             extent[axis]=minmax[axis][1]-minmax[axis][0]
             center[axis]=extent[axis]/2.0+minmax[axis][0]
         return [center['x'],center['y'],center['z']]
@@ -1852,9 +1852,9 @@
             if found==len(atomlist):
                 break
         if abs(totphi)<0.01 or abs(totcrg)<0.01:
-            print 'total abs phi',totphi
-            print 'total abs crg',totcrg
-            print 'net charge   ',netcrg
+            print('total abs phi',totphi)
+            print('total abs crg',totcrg)
+            print('net charge   ',netcrg)
             PDB2PKAError( 'Something is rotten')
 
         return energy
@@ -1881,7 +1881,7 @@
                 # Yes!
                 #
                 with open(result_file,'rb') as fd:
-                    potentials=cPickle.load(fd)
+                    potentials=pickle.load(fd)
                     loaded=True
         #
         # Run calc again if needed
@@ -1895,7 +1895,7 @@
                 self.APBS=None
             if save_results:
                 with open(result_file,'wb') as fd:
-                    cPickle.dump(potentials,fd)
+                    pickle.dump(potentials,fd)
 
         return potentials
 
@@ -1932,7 +1932,7 @@
                 text = "Could not find radius for atom %s" % atomname
                 text += " in residue %s %i" % (residue.name, residue.resSeq)
                 text += " while attempting to set radius!"
-                raise ValueError, text
+                raise ValueError(text)
     #
     # ------------------------------------
     #
@@ -1959,7 +1959,7 @@
                 text = "Could not find charge for atom %s" % atomname
                 text += " in residue %s %i" % (residue.name, residue.resSeq)
                 text += " while attempting to set charge!"
-                raise ValueError, text
+                raise ValueError(text)
         return
     #
     # ----------------------------
@@ -2038,8 +2038,8 @@
             charge, radius = self.forcefield.getParams1(residue, atomname)
             initialmap[atomname] = charge
             if charge is None:
-                print atomname,charge
-                print residue.isCterm
+                print(atomname,charge)
+                print(residue.isCterm)
                 raise PDB2PKAError('Charge on atom is None')
             sum+=charge
         if abs(sum)<0.001:
@@ -2061,7 +2061,7 @@
 
                 charge, radius = self.forcefield.getParams1(residue, atomname)
                 sum=sum+charge
-                if initialmap.has_key(atomname):
+                if atomname in initialmap:
                     initcharge = initialmap[atomname]
                     if charge != initcharge:
                         if not atomname in atomnames:
@@ -2072,7 +2072,7 @@
             #
             # Check that no atoms were removed
             #
-            for atom in initialmap.keys():
+            for atom in list(initialmap.keys()):
                 if not atom in residue.get('map'):
                     atomnames.append(atom)
         #
@@ -2148,8 +2148,8 @@
             # Did we add anything?
             #
             if added is None and sum>0.001:
-                print sum
-                print atomnames
+                print(sum)
+                print(atomnames)
                 PDB2PKAError('Could not find integer charge state')
         #
         # Did we just want a neutral state identification?
@@ -2162,7 +2162,7 @@
         # No, we wanted the atomnames
         #
         if atomnames==[]:
-            print 'Did not find any atoms for ',residue.resSeq
+            print('Did not find any atoms for ',residue.resSeq)
             PDB2PKAError('Something wrong with charges')
 
         for atomname in atomnames:
@@ -2192,7 +2192,7 @@
         self.routines.write("Finding Titratable groups....\n")
         sys.stdout.flush()
         #
-        pKagroupList=self.pKagroups.keys()
+        pKagroupList=list(self.pKagroups.keys())
         #
         for chain in self.protein.getChains():
             for residue in chain.get("residues"):
@@ -2263,7 +2263,7 @@
         if amb == None:
             text = "Could not find hydrogen ambiguity "
             text += "for titratable group %s!" % group
-            raise ValueError, text
+            raise ValueError(text)
         return amb
 
     #
@@ -2288,7 +2288,7 @@
                          'CTR01c': '1', 'CTR01t': '2', 'CTR02c': '3', 'CTR02t': '4', 'CTR-': '0'}
         filename = TITRATIONFILE
         if not os.path.isfile(TITRATIONFILE):
-            raise ValueError, "Could not find TITRATION.DAT!"
+            raise ValueError("Could not find TITRATION.DAT!")
 
         titration_file = open(filename)
 
@@ -2306,18 +2306,18 @@
                 line = titration_file.readline()
                 if line[:8] != 'Residue:':
                     text = "Wrong line found when looking for 'Residue'"
-                    raise ValueError, "%s: %s" % (text, line)
+                    raise ValueError("%s: %s" % (text, line))
 
                 resname = string.strip(string.split(line)[1])
 
                 line = titration_file.readline()
                 if line[:10] != 'Grouptype:':
                     text = "Wrong line found when looking for 'Grouptype'"
-                    raise ValueError, "%s: %s" % (text, line)
+                    raise ValueError("%s: %s" % (text, line))
 
                 type = string.lower(string.strip(string.split(line)[1]))
                 if type != 'acid' and type != 'base':
-                    raise ValueError, 'Group type must be acid or base!'
+                    raise ValueError('Group type must be acid or base!')
 
                 line = titration_file.readline()
                 while 1:
@@ -2334,7 +2334,7 @@
 
                     if line[:11] != 'Transition:':
                         text = "Wrong line found when looking for 'Transition:'"
-                        raise ValueError, "%s: %s" % (text, line)
+                        raise ValueError("%s: %s" % (text, line))
 
                     split=string.split(line[11:],'->')
                     for number in string.split(split[0], ','):
@@ -2353,7 +2353,7 @@
                     #
                     if line[:10]!='Model_pKa:':
                         text = "Wrong line found when looking for 'Model_pKa'"
-                        raise ValueError, "%s: %s" % (text, line)
+                        raise ValueError("%s: %s" % (text, line))
 
                     modelpKa = float(string.split(line)[1])
 
@@ -2401,7 +2401,7 @@
 #
 
 def smooth(xdiel,ydiel,zdiel):
-    print '\nSmooting dielectric constant using Gaussian filter:\n'
+    print('\nSmooting dielectric constant using Gaussian filter:\n')
 
     diel=[xdiel,ydiel,zdiel]
     for d in diel:
@@ -2441,18 +2441,18 @@
 
 
 
-    print "These should pass without issue"
-    print "Run acid curve"
+    print("These should pass without issue")
+    print("Run acid curve")
     routines = pKaRoutines(None, None, None, None, '', maps = None, sd =None,
                  restart=False, pairene=1.0, test_mode=True)
     routines.find_pH_at_0_5('zero_to_neg_one curve base', zero_to_neg_one, False)
 
-    print "Run base curve"
+    print("Run base curve")
     routines.find_pH_at_0_5('one_to_zero curve acid', one_to_zero, True)
 
-    print "These should print warnings"
+    print("These should print warnings")
     all_zero_curve = dict((ph,0.0) for ph in ph_list)
-    print "Run all zero curves"
+    print("Run all zero curves")
     routines.find_pH_at_0_5('All zero curve acid', all_zero_curve, False)
     routines.find_pH_at_0_5('All zero curve base', all_zero_curve, True)
 
@@ -2492,9 +2492,9 @@
 
     routines.find_pH_at_0_5('negative_interpolation_curve acid', negative_interpolation_curve, False)
 
-    print 'Accumulated warnings:'
+    print('Accumulated warnings:')
     pprint(routines.warnings)
-    print 'ph values:'
+    print('ph values:')
     pprint(routines.ph_at_0_5)
 
 
--- a/pdb2pka/prepare_pKa_ligand.py
+++ b/pdb2pka/prepare_pKa_ligand.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
 
 #
 # Take a mol2 file as input, convert to SMILES, search pka_lig_tool database and produce
@@ -70,8 +70,8 @@
         
     def read_mol2(self,mol2lines):
         """Parse the mol2 lines"""
-        import StringIO, string
-        mol2fileobj=StringIO.StringIO(string.join(mol2lines))
+        import io, string
+        mol2fileobj=io.StringIO(string.join(mol2lines))
         #
         # Use the mol2 parser in pdb2pqr
         #
@@ -121,18 +121,18 @@
 
     def search_pka_ligtool(self,smiles):
         """Search the pka_lig_tool database for a ligand match"""
-        import StringIO, urllib
+        import io, urllib.request, urllib.parse, urllib.error
     
         #
         # Get the XML data from the server
         #
-        args=urllib.urlencode({'smiles':smiles})
+        args=urllib.parse.urlencode({'smiles':smiles})
         thisurl= url %(server,args)
-        f=urllib.urlopen(thisurl)
+        f=urllib.request.urlopen(thisurl)
         text=f.read()
-        output = StringIO.StringIO(text)
-        print text
-        print smiles
+        output = io.StringIO(text)
+        print(text)
+        print(smiles)
         return text
         #
         # Parse the XML
@@ -141,39 +141,39 @@
         xmldoc = minidom.parse(output)
 
         ligands = xmldoc.firstChild
-        print 'Search-type was: %s'%ligands.attributes['Type'].value
+        print('Search-type was: %s'%ligands.attributes['Type'].value)
         for ligand in ligands.childNodes:
             atoms = ligand.getElementsByTagName('Atoms')[0]
             mol2 = ligand.getElementsByTagName('mol2')[0]
 
-            print '*'*50
-            print 'Found ligand \"%s\"'%ligand.attributes['Name'].value
-            print
-            print 'mol2 file'
-            print '-'*50
-            print mol2.firstChild.data
-            print
-            print ' Atoms and associated pKa values'
-            print '-'*50
+            print('*'*50)
+            print('Found ligand \"%s\"'%ligand.attributes['Name'].value)
+            print()
+            print('mol2 file')
+            print('-'*50)
+            print(mol2.firstChild.data)
+            print()
+            print(' Atoms and associated pKa values')
+            print('-'*50)
             for atom in atoms.childNodes:
-                print '%4s %4s %-6s'%(atom.attributes['Name'].value,
+                print('%4s %4s %-6s'%(atom.attributes['Name'].value,
                                      atom.attributes['Number'].value,
-                                     atom.attributes['Type'].value)
+                                     atom.attributes['Type'].value))
                 pkas = atom.firstChild
                 if pkas:
                     for pka in pkas.childNodes:
-                        print '          Value:            ',pka.attributes['Value'].value
-                        print '          Temperature:      ',pka.attributes['Temperature'].value
-                        print '          pH:               ',pka.attributes['pH'].value
-                        print '          Solvent:          ',pka.attributes['Solvent'].value
-                        print '          Salt type:        ',pka.attributes['Salt_type'].value
-                        print '          Salt conc.:       ',pka.attributes['Salt_conc.'].value
-                        print '          Titratable group: ',pka.attributes['Titratable_group'].value
-                        print '          Most bio. rel.:   ',pka.attributes['Most_bio._relavent'].value
-                        print '          Reference:        ',pka.attributes['Reference'].value
-                        print '          Comment:          ',pka.attributes['Comment'].value
+                        print('          Value:            ',pka.attributes['Value'].value)
+                        print('          Temperature:      ',pka.attributes['Temperature'].value)
+                        print('          pH:               ',pka.attributes['pH'].value)
+                        print('          Solvent:          ',pka.attributes['Solvent'].value)
+                        print('          Salt type:        ',pka.attributes['Salt_type'].value)
+                        print('          Salt conc.:       ',pka.attributes['Salt_conc.'].value)
+                        print('          Titratable group: ',pka.attributes['Titratable_group'].value)
+                        print('          Most bio. rel.:   ',pka.attributes['Most_bio._relavent'].value)
+                        print('          Reference:        ',pka.attributes['Reference'].value)
+                        print('          Comment:          ',pka.attributes['Comment'].value)
 
-            print '*'*50
+            print('*'*50)
         return
 
     def get_allhyd_state(self):
@@ -187,10 +187,10 @@
 
 
 if __name__=='__main__':
-    print
-    print 'Get pKa values and structures of protonation states for a ligand'
-    print 'Chresten Soendergaard, Paul Czodrowski, Jens Erik Nielsen 2006-2010'
-    print
+    print()
+    print('Get pKa values and structures of protonation states for a ligand')
+    print('Chresten Soendergaard, Paul Czodrowski, Jens Erik Nielsen 2006-2010')
+    print()
     import sys, os
     from optparse import OptionParser
     parser = OptionParser(usage='%prog [options] <file>',version='%prog 1.0')
@@ -228,7 +228,7 @@
                         break
                 #
                 if mol2file:
-                    print mol2file
+                    print(mol2file)
                     fd=open(mol2file)
                     mol2lines=fd.readlines()
                     fd.close()
@@ -241,11 +241,11 @@
                     except:
                         import sys
                         failed.append([mol2file,sys.exc_info()[0]])
-                        print 'FAILED'
-                        print sys.exc_info()[0]
-        print failed
-        print 'OK',len(ok)
-        print 'FAILED',len(failed)
+                        print('FAILED')
+                        print(sys.exc_info()[0])
+        print(failed)
+        print('OK',len(ok))
+        print('FAILED',len(failed))
                         
 
     
--- a/pka.py
+++ b/pka.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
 #
 # pKa calculations with APBS
 #
@@ -28,11 +28,11 @@
     #
     # Print usage guidelines
     #
-    print 'Usage: pka.py --ff <forcefield> --lig <ligand in MOL2> --pdie <protein diel cons> --maps <1 for using provided 3D maps; 2 for genereting new maps>'
-    print '--xdiel <xdiel maps> --ydiel <ydiel maps> --zdiel <zdiel maps> --kappa <ion-accessibility map> '
-    print '--smooth <st.dev [A] of Gaussian smooting of 3D maps at the boundary, bandthwith=3 st.dev> <pdbfile>'
-    print 'Force field can be amber, charmm and parse'
-    print
+    print('Usage: pka.py --ff <forcefield> --lig <ligand in MOL2> --pdie <protein diel cons> --maps <1 for using provided 3D maps; 2 for genereting new maps>')
+    print('--xdiel <xdiel maps> --ydiel <ydiel maps> --zdiel <zdiel maps> --kappa <ion-accessibility map> ')
+    print('--smooth <st.dev [A] of Gaussian smooting of 3D maps at the boundary, bandthwith=3 st.dev> <pdbfile>')
+    print('Force field can be amber, charmm and parse')
+    print()
     return
 
 #
@@ -48,9 +48,9 @@
             routines:   The routines object as generated by PDB2PQR
             forcefield: The forcefield object as generated by PDB2PQR
     """
-    print
-    print 'PDB2PQR pKa calculations'
-    print
+    print()
+    print('PDB2PQR pKa calculations')
+    print()
 
     parser = optparse.OptionParser()
 
@@ -184,7 +184,7 @@
         try:
             ligand = open(options.ligand, 'rU')
         except IOError:
-            print 'Unable to find ligand file %s! Skipping...' % options.ligand
+            print('Unable to find ligand file %s! Skipping...' % options.ligand)
 
     #Set up the protien object
     #In the standalone version of pdb2pka this is redundent but needed so we emulate the
@@ -193,8 +193,8 @@
     pdbfile = getPDBFile(input_path)
     pdblist, errlist = readPDB(pdbfile)
     if len(errlist) != 0 and verbose:
-        print "Warning: %s is a non-standard PDB file.\n" %input_path
-        print errlist
+        print("Warning: %s is a non-standard PDB file.\n" %input_path)
+        print(errlist)
     #
     # Read the definition file
     #
@@ -278,7 +278,7 @@
     pdblist, errlist = readPDB(pdbfile)
 
     if verbose:
-        print "Beginning PDB2PKA...\n"
+        print("Beginning PDB2PKA...\n")
     #
     # Read the definition file
     #
@@ -304,9 +304,9 @@
     # Print something for some reason?
     #
     if verbose:
-        print "Created protein object -"
-        print "\tNumber of residues in protein: %s" % myProtein.numResidues()
-        print "\tNumber of atoms in protein   : %s" % myProtein.numAtoms()
+        print("Created protein object -")
+        print("\tNumber of residues in protein: %s" % myProtein.numResidues())
+        print("\tNumber of atoms in protein   : %s" % myProtein.numAtoms())
     #
     # Set up all other routines
     #
@@ -357,8 +357,8 @@
                 templist = []
                 Lig.make_up2date(residue)
                 net_charge=0.0
-                print 'Ligand',residue
-                print 'Atom\tCharge\tRadius'
+                print('Ligand',residue)
+                print('Atom\tCharge\tRadius')
                 for atom in residue.getAtoms():
                     if atom.mol2charge:
                         atom.ffcharge=atom.mol2charge
@@ -372,7 +372,7 @@
                     # Assign radius
                     #
                     atom.radius = Lig.ligand_props[atom.name]["radius"]
-                    print '%s\t%6.4f\t%6.4f' %(atom.name,atom.ffcharge,atom.radius)
+                    print('%s\t%6.4f\t%6.4f' %(atom.name,atom.ffcharge,atom.radius))
                     if atom in misslist:
                         misslist.pop(misslist.index(atom))
                         templist.append(atom)
@@ -397,7 +397,7 @@
                 #
                 # Print the net charge
                 #
-                print 'Net charge for ligand %s is: %5.3f' %(residue.name,net_charge)
+                print('Net charge for ligand %s is: %5.3f' %(residue.name,net_charge))
         #
         # Temporary fix; if ligand was successful, pull all ligands from misslist
         # Not sure if this is needed at all here ...? (Jens wrote this)
@@ -410,9 +410,9 @@
                 misslist.remove(atom)
 
     if verbose:
-        print "Created protein object (after processing myRoutines) -"
-        print "\tNumber of residues in protein: %s" % myProtein.numResidues()
-        print "\tNumber of atoms in protein   : %s" % myProtein.numAtoms()
+        print("Created protein object (after processing myRoutines) -")
+        print("\tNumber of residues in protein: %s" % myProtein.numResidues())
+        print("\tNumber of atoms in protein   : %s" % myProtein.numAtoms())
     #
     # Create the APBS input file
     #
@@ -427,11 +427,11 @@
     # For convenience
     #
     igen.pdie = pdie
-    print 'Setting protein dielectric constant to ',igen.pdie
+    print('Setting protein dielectric constant to ',igen.pdie)
     igen.sdie=sdie
     igen.maps=maps
     if maps==1:
-        print "Using dielectric and mobile ion-accessibility function maps in PBE"
+        print("Using dielectric and mobile ion-accessibility function maps in PBE")
         if xdiel:
             igen.xdiel = xdiel
         else:
@@ -445,14 +445,14 @@
         else:
             raise PDB2PKAError("Z dielectric map is missing\n")
 
-        print 'Setting dielectric function maps: %s, %s, %s'%(igen.xdiel,igen.ydiel,igen.zdiel)
+        print('Setting dielectric function maps: %s, %s, %s'%(igen.xdiel,igen.ydiel,igen.zdiel))
 
         if kappa:
             igen.kappa = kappa
         else:
             raise PDB2PKAError("Mobile ion-accessibility map is missing\n")
 
-        print 'Setting mobile ion-accessibility function map to: ',igen.kappa
+        print('Setting mobile ion-accessibility function map to: ',igen.kappa)
 
         if sd:
             xdiel_smooth, ydiel_smooth, zdiel_smooth = smooth(xdiel,ydiel,zdiel)
@@ -475,6 +475,6 @@
     mypkaRoutines = pka_routines.pKaRoutines(protein, routines, forcefield, apbs_setup, output_dir, maps, sd,
                                              restart=not options.resume, pairene=options.pairene)
 
-    print 'Doing full pKa calculation'
+    print('Doing full pKa calculation')
     mypkaRoutines.runpKa()
 
--- a/propka30/Source/bonds.py
+++ b/propka30/Source/bonds.py
@@ -38,7 +38,7 @@
 #-------------------------------------------------------------------------------------------------------
 import pickle,sys,os,math
 
-from lib import pka_print
+from .lib import pka_print
 
 
 class bondmaker:
@@ -250,9 +250,9 @@
         
         # apply table information on pi-electrons
         for atom in ligand.atoms:
-            if atom.name in self.number_of_pi_electrons_in_bonds_ligands.keys():
+            if atom.name in list(self.number_of_pi_electrons_in_bonds_ligands.keys()):
                 atom.number_of_pi_electrons_in_double_and_triple_bonds = self.number_of_pi_electrons_in_bonds_ligands[atom.name]
-            if atom.name in self.number_of_pi_electrons_in_conjugate_bonds_in_ligands.keys():
+            if atom.name in list(self.number_of_pi_electrons_in_conjugate_bonds_in_ligands.keys()):
                 atom.number_of_pi_electrons_in_conjugate_double_and_triple_bonds = self.number_of_pi_electrons_in_conjugate_bonds_in_ligands[atom.name]
 
             # just in case any protein residues are included
@@ -360,7 +360,7 @@
             
 
         # assign bonds 
-        keys = self.boxes.keys()
+        keys = list(self.boxes.keys())
         for key in keys:
             self.find_bonds_for_atoms(self.boxes[key])
             
@@ -374,7 +374,7 @@
             for by in [y-1,y,y+1]:
                 for bz in [z-1,z,z+1]:   
                     key = self.box_key(bx,by,bz)
-                    if key in self.boxes.keys():
+                    if key in list(self.boxes.keys()):
                         self.boxes[key].append(atom)
 
                         #pka_print(atom,'->',key,':',len(self.boxes[key]))
--- a/propka30/Source/calculator.py
+++ b/propka30/Source/calculator.py
@@ -38,7 +38,7 @@
 #-------------------------------------------------------------------------------------------------------
 import math, random, string
 
-from lib import pka_print
+from .lib import pka_print
 
 
 def InterAtomDistance(atom1, atom2):
@@ -218,7 +218,7 @@
     dV            = 0.00
     volume        = 0.00
     min_distance_4th = pow(2.75, 4)
-    for chainID in atoms.keys():
+    for chainID in list(atoms.keys()):
       for key in atoms[chainID]["keys"]:
         for atom in atoms[chainID][key]:
           if atom.element != "H":
@@ -273,7 +273,7 @@
       local_cutoff = 0.00
     residue.Nmass = 0
     residue.Nlocl = 0
-    for chainID in atoms.keys():
+    for chainID in list(atoms.keys()):
       for key in atoms[chainID]["keys"]:
         for atom in atoms[chainID][key]:
           if atom.element != "H":
@@ -312,7 +312,7 @@
     Nlocl_his6   = 0
     residue.Nmass = 0
     residue.Nlocl = 0
-    for chainID in atoms.keys():
+    for chainID in list(atoms.keys()):
       for key in atoms[chainID]["keys"]:
         for atom in atoms[chainID][key]:
           HYDROGEN_ATOM = ( (atom.name[0] == 'H') or (atom.name[0] in string.digits and atom.name[1] == 'H') )
--- a/propka30/Source/chain.py
+++ b/propka30/Source/chain.py
@@ -37,10 +37,10 @@
 #   Journal of Chemical Theory and Computation, 7, 525-537 (2011)
 #-------------------------------------------------------------------------------------------------------
 import math, sys
-import lib
+from . import lib
 pka_print = lib.pka_print
-import mutate
-from residue import Residue
+from . import mutate
+from .residue import Residue
 
 
 class Chain:
--- a/propka30/Source/compare.py
+++ b/propka30/Source/compare.py
@@ -37,7 +37,7 @@
 #   Journal of Chemical Theory and Computation, 7, 525-537 (2011)
 #-------------------------------------------------------------------------------------------------------
 import  sys, math, string
-from lib import int2roman, convertResidueCode, pka_print
+from .lib import int2roman, convertResidueCode, pka_print
 
 
 def compareFoldingContributions(target=None, template=None, options=None):
@@ -48,7 +48,7 @@
     4. calculate the difference
     5. printout sorted result
     """
-    from mutate import readAlignmentFiles
+    from .mutate import readAlignmentFiles
 
     # checking that pKa values are available
     checkDonePKA(target, template)
@@ -77,7 +77,7 @@
     make a list of pdbcodes in 'alignment'
     """
     names = [name]
-    for key in alignment.keys():
+    for key in list(alignment.keys()):
       if key not in names:
         names.append(key)
 
@@ -115,8 +115,8 @@
 """
     str = str[:-1]
     pka_print(str)
-    for key1 in alignment.keys():
-      for key2 in alignment[key1].keys():
+    for key1 in list(alignment.keys()):
+      for key2 in list(alignment[key1].keys()):
         str = "    %s 100 %s 0" % (key2, "%")
         index = 0
         for code in alignment[key1][key2]['sequence']:
--- a/propka30/Source/corresponding_atoms.py
+++ b/propka30/Source/corresponding_atoms.py
@@ -40,7 +40,7 @@
 #-------------------------------------------------------------------------------------------------------
 
 import sys
-from lib import residueList, atomList, pka_print
+from .lib import residueList, atomList, pka_print
 
 
 #    NOTE:
--- a/propka30/Source/coupled_residues.py
+++ b/propka30/Source/coupled_residues.py
@@ -38,7 +38,7 @@
 #-------------------------------------------------------------------------------------------------------
 
 import math
-from lib import pka_print
+from .lib import pka_print
 
 max_intrinsic_pKa_diff = 2.0
 min_interaction_energy = 0.5
--- a/propka30/Source/debug.py
+++ b/propka30/Source/debug.py
@@ -37,10 +37,10 @@
 #   Journal of Chemical Theory and Computation, 7, 525-537 (2011)
 #-------------------------------------------------------------------------------------------------------
 import string
-import lib
-import calculator as calculate
+from . import lib
+from . import calculator as calculate
 
-from lib import pka_print
+from .lib import pka_print
 
 
 def interactionMatrix(interaction):
@@ -69,9 +69,9 @@
     printing out all information in resInfo
     """
     pka_print("in resInfo:")
-    for key1 in resInfo.keys():
+    for key1 in list(resInfo.keys()):
       pka_print(" --- %s ---" % (key1))
-      for key2 in resInfo[key1].keys():
+      for key2 in list(resInfo[key1].keys()):
         pka_print(key2, resInfo[key1][key2])
 
 
@@ -162,9 +162,9 @@
     """
     Prints out alignment information for debugging
     """
-    for key in alignment.keys():
+    for key in list(alignment.keys()):
       pka_print( " --- %s ---" % (key) )
-      for key2 in alignment[key].keys():
+      for key2 in list(alignment[key].keys()):
         pka_print("%s %5d%2s" % (alignment[key][key2]["name"], alignment[key][key2]["resNumb"], alignment[key][key2]["chainID"]))
         pka_print("%s\n" % (alignment[key][key2]["sequence"]))
 
--- a/propka30/Source/determinants.py
+++ b/propka30/Source/determinants.py
@@ -38,12 +38,12 @@
 #-------------------------------------------------------------------------------------------------------
 import math, time
 
-import iterative
-import lib
-from lib import pka_print
+from . import iterative
+from . import lib
+from .lib import pka_print
 #import debug
-import calculator as calculate
-from   determinant import Determinant
+from . import calculator as calculate
+from   .determinant import Determinant
 
 
 def setDeterminants(propka_residues, version=None, options=None):
--- a/propka30/Source/iterative.py
+++ b/propka30/Source/iterative.py
@@ -37,11 +37,11 @@
 #   Journal of Chemical Theory and Computation, 7, 525-537 (2011)
 #-------------------------------------------------------------------------------------------------------
 import math, time
-import calculator as calculate
-import lib
+from . import calculator as calculate
+from . import lib
 pka_print = lib.pka_print
 #import debug
-from determinant import Determinant
+from .determinant import Determinant
 
 
 # Some library functions for the interative pKa determinants
--- a/propka30/Source/lib.py
+++ b/propka30/Source/lib.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 #
 # * This library is free software; you can redistribute it and/or
 # * modify it under the terms of the GNU Lesser General Public
@@ -40,7 +40,7 @@
 
 
 import string, sys, copy, math
-import output
+from . import output
 
 
 def loadOptions():
@@ -199,7 +199,7 @@
           if filename not in options.alignment:
             options.alignment.append( filename )
         else:
-          for code in mutation.keys():
+          for code in list(mutation.keys()):
             filename = "%s.pir" % ( extractName(code) )
             if filename not in options.alignment:
               options.alignment.append( filename )
@@ -762,7 +762,7 @@
     numerals = { 1 : "I", 4 : "IV", 5 : "V", 9 : "IX", 10 : "X", 40 : "XL",
         50 : "L", 90 : "XC", 100 : "C", 400 : "CD", 500 : "D", 900 : "CM", 1000 : "M" }
     result = ""
-    for value, numeral in sorted(numerals.items(), reverse=True):
+    for value, numeral in sorted(list(numerals.items()), reverse=True):
         while number >= value:
             result += numeral
             number -= value
--- a/propka30/Source/ligand.py
+++ b/propka30/Source/ligand.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 #
 # * This library is free software; you can redistribute it and/or
 # * modify it under the terms of the GNU Lesser General Public
@@ -40,7 +40,7 @@
 
 
 import sys, pdb, protonate, lib, bonds
-from vector_algebra import *
+from .vector_algebra import *
 pka_print = lib.pka_print
 
 all_sybyl_types = [
@@ -185,7 +185,7 @@
             atom.residue = self
 
         #self.remove_ions()
-        self.configuration_keys = atoms[0].configurations.keys()
+        self.configuration_keys = list(atoms[0].configurations.keys())
         
         # create ligand residue objects
         self.ligand_residues = []
@@ -265,7 +265,7 @@
         for atom in self.atoms:
             # check if we already have assigned a name to this atom
             if hasattr(atom, 'sybyl_assigned'):
-                print(atom.resName, atom.numb, atom.name, 'alreadyassigned')
+                print((atom.resName, atom.numb, atom.name, 'alreadyassigned'))
                 continue
 
             # find some properties of the atom
@@ -312,7 +312,7 @@
 
             if atom.get_element() == 'C':
                 # check for amide
-                if 'O' in bonded_elements.values() and 'N' in bonded_elements.values():
+                if 'O' in list(bonded_elements.values()) and 'N' in list(bonded_elements.values()):
                     self.set_type(atom, 'C.2')
                     for b in atom.bonded_atoms:
                         if b.get_element() == 'N':
--- a/propka30/Source/mutate.py
+++ b/propka30/Source/mutate.py
@@ -38,10 +38,10 @@
 #-------------------------------------------------------------------------------------------------------
 import math, os, sys, re
 
-import lib
+from . import lib
 pka_print = lib.pka_print
-import output
-import pdb
+from . import output
+from . import pdb
 #import debug
 
 
@@ -119,7 +119,7 @@
     best_mutation = "WT"
     dG_ref = proteins["WT"][0]
     pka_print("%8s %6s\n%8.2lf %6.2lf  %s" % ("dG_fold", "ddG_fold", dG_ref, dG_ref-dG_ref, "WT"))
-    for label in proteins.keys():
+    for label in list(proteins.keys()):
       if label != "WT":
         dG_mut, protein, full_mutation = proteins[label]
         pka_print("%8.2lf %6.2lf  %s" % (dG_mut, dG_mut-dG_ref, label))
@@ -163,7 +163,7 @@
     """
     returning a new protein, mutation approach determined from options.mutator
     """
-    from protein import Protein
+    from .protein import Protein
     pka_print("mutator: %s" % (options.mutator.label))
     pka_print("type: %s\n" % (options.mutator.type))
 
@@ -174,9 +174,9 @@
         alignment = readAlignmentFiles(filenames=options.alignment, mesophile=protein.name, options=options)
       pdbcode, rest = splitStringMutationInTwo(mutation)
       if pdbcode == None:
-        mutation = remakeStringMutation(mutation, protein.atoms.keys())
+        mutation = remakeStringMutation(mutation, list(protein.atoms.keys()))
       else:
-        mutation = remakeStringMutation(mutation, protein.atoms.keys(), atoms[pdbcode].keys())
+        mutation = remakeStringMutation(mutation, list(protein.atoms.keys()), list(atoms[pdbcode].keys()))
       dict_mutation = recastIntoDictionary(mutation, alignment=alignment, protein=protein)
     elif isinstance(mutation, dict):
       # this is already in dictionary-format
@@ -233,7 +233,7 @@
     print out the mutation in formatted form
     """
     pka_print(" ----- mutation ----- \n", "  residue    rmsd")
-    for key in mutation.keys():
+    for key in list(mutation.keys()):
       pka_print("  %s %6.2lf" % (mutation[key]['label'], mutation[key]['rmsd']))
 
 
@@ -325,7 +325,7 @@
     
     
         # store it as alignment-pairs in 'alignments'
-        for key in alignment.keys():
+        for key in list(alignment.keys()):
           if re.search(mesophile, key):
             """ do nothing, this matches the mesophile """
             # pka_print("match,  don't do \"%s\" \"%s\":" % (mesophile, key))
@@ -359,7 +359,7 @@
     get the translation for the 'back-tracking' in the 'alignment-mutation'
     """
     rmsd = 0.00; number_of_atoms = 0; translation = [0.00, 0.00, 0.00]
-    for key in backbone.keys():
+    for key in list(backbone.keys()):
       number_of_atoms += 1
       atom1, atom2 = backbone[key]
       dX = atom2.x - atom1.x;
@@ -381,7 +381,7 @@
     creating a copy of the atoms dictionary
     """
     newAtoms = {}
-    for chainID in atoms.keys():
+    for chainID in list(atoms.keys()):
       # creating new chain dictionary
       newAtoms[chainID] = {'keys': []}
       for key in atoms[chainID]['keys']:
@@ -402,7 +402,7 @@
     newAtoms = copyAtomsDictionary(atoms=protein.atoms, options=options)
     configs  = [ protein.configurations[0] ]
 
-    for key in mutation.keys():
+    for key in list(mutation.keys()):
 
       singleMutateAtomsDictionary(atoms=newAtoms, 
                                   template=atoms, 
@@ -476,7 +476,7 @@
       newAtoms[chainID] = {"keys": []}
       for key in protein.atoms[chainID]["keys"]:
         atomList = []
-        if key in mutation.keys():
+        if key in list(mutation.keys()):
           # adding residue from thermophile data
           code = mutation[key]['pdb']
           key1 = mutation[key]['template']
@@ -530,7 +530,7 @@
     elif isinstance(generic_mutation, list):
       mutation_list = generic_mutation
     elif isinstance(generic_mutation, dict):
-      for key in generic_mutation.keys():
+      for key in list(generic_mutation.keys()):
         chainID1 = generic_mutation[key]['target'][-1]
         code1, resName = lib.convertResidueCode(resName=generic_mutation[key]['target'][:3])
         resNumb = int(generic_mutation[key]['target'][3:7])
@@ -718,7 +718,7 @@
           if   isinstance(mutations, dict):
             mutation[ keys[i] ] = mutations[ keys[i] ]
           elif isinstance(mutations, list):
-            for key in mutations[i].keys():
+            for key in list(mutations[i].keys()):
               mutation[ key ] = mutations[i][ key ]
         else:
           short += "N"
@@ -826,7 +826,7 @@
 
     label = ""
     if isinstance(generic_mutation, dict):
-      for key in generic_mutation.keys():
+      for key in list(generic_mutation.keys()):
         code1, resName = lib.convertResidueCode(resName=key[:3])
         code2, resName = lib.convertResidueCode(resName=generic_mutation[key]['label'][:3])
         label += "/%s%d%s" % (code1, int(key[3:7]), code2)
@@ -900,7 +900,7 @@
         new[target_label]['target']   = target_label
       mutation = new
     elif isinstance(mutation, dict):
-      for key in mutation.keys():
+      for key in list(mutation.keys()):
         mutation[key]['target'] = key
         if "label" not in mutation[key] and "template" in mutation[key]:
           # generating new label
@@ -917,7 +917,7 @@
     rmsd_new = 0.00
     rmsd_old = 0.00
     number_of_points = len(position['target'])
-    for key in position['target'].keys():
+    for key in list(position['target'].keys()):
       for i in range(3):
         rmsd_new += pow( (position['new'][key][i] - position['target'][key][i]), 2)
         rmsd_old += pow( (position['old'][key][i] - position['target'][key][i]), 2)
@@ -1010,8 +1010,8 @@
     """
     This routine overlaps two residues based on an array of atom labels, 'center'
     """
-    from corresponding_atoms import makeCorrespondingAtomNames
-    from rotate import rotatePosition, translatePosition, makeCrossProduct, calculateVectorLength, \
+    from .corresponding_atoms import makeCorrespondingAtomNames
+    from .rotate import rotatePosition, translatePosition, makeCrossProduct, calculateVectorLength, \
                               makeScalarProduct, generateRandomDisplacement, generateRandomRotation, rotateAtoms, translateAtoms
 
     # create copy of target
@@ -1022,7 +1022,7 @@
     dR       = [None, None, None]
     center   = ['CA']
 
-    for key in mutation.keys():
+    for key in list(mutation.keys()):
 
       # mutate the new atoms dictionary (position unchanged)
       singleMutateAtomsDictionary(atoms=newAtoms,
@@ -1076,7 +1076,7 @@
       # initializing iterations: copy over 'template' to 'new' & 'old'
       chainNumb = ord("A")
       center    = []
-      for name in position['template'].keys():
+      for name in list(position['template'].keys()):
         center.append(name)
         position['new'][name] = [position['template'][name][0], position['template'][name][1], position['template'][name][2]]
         position['old'][name] = [position['template'][name][0], position['template'][name][1], position['template'][name][2]]
@@ -1101,7 +1101,7 @@
         if rmsd_new < rmsd_old:
           # pka_print("REMARK   update structure (iter %3d (%8.3lf%8.3lf))" % (iter, rmsd_new, rmsd_old))
           chainNumb += 1
-          for name in position['old'].keys():
+          for name in list(position['old'].keys()):
             for i in range(3):
               position['old'][name][i] = position['new'][name][i]
             if False:
@@ -1111,7 +1111,7 @@
           rotateAtoms(newAtoms[target_chainID][label], axis, theta, center=center)
           translateAtoms(newAtoms[target_chainID][label], dR)
         else:
-          for name in position['new'].keys():
+          for name in list(position['new'].keys()):
             for i in range(3):
               position['new'][name][i] = position['old'][name][i]
 
@@ -1172,7 +1172,7 @@
               atoms[pdbcode] = pdb.readPDB(filename=pdbcode)
           else:
             # extracting pdbcode if mutation is in dictionary format
-            for key in mutation.keys():
+            for key in list(mutation.keys()):
               pdbcode = mutation[key]['pdb']
               if pdbcode not in atoms:
                 atoms[pdbcode] = pdb.readPDB(filename=pdbcode)
@@ -1201,10 +1201,10 @@
     translate = dictionary for translating target and template residues to the origin
     center    = array with atom names in 'position-dictionary': used for 'rotation-center'
     """
-    from rotate import rotatePosition, translatePosition, makeCrossProduct, calculateVectorLength, makeScalarProduct, generateRandomDisplacement, generateRandomRotation
+    from .rotate import rotatePosition, translatePosition, makeCrossProduct, calculateVectorLength, makeScalarProduct, generateRandomDisplacement, generateRandomRotation
 
     # start mutations
-    for label in mutation.keys():
+    for label in list(mutation.keys()):
 
       pka_print("%s ==> %s" % (label, mutation[label]['template']))
 
@@ -1219,10 +1219,10 @@
 
       # debug printout
       if False:
-        for key in position['target'].keys():
+        for key in list(position['target'].keys()):
           Rx, Ry, Rz = position['target'][key]
           pka_print("ATOM    669  %-3s %s A  89    %8.3lf%8.3lf%8.3lf" % (key, original_residue.resName, Rx, Ry, Rz))
-        for key in position['template'].keys():
+        for key in list(position['template'].keys()):
           Rx, Ry, Rz = position['template'][key]
           pka_print("ATOM    669  %-3s GLU B  89    %8.3lf%8.3lf%8.3lf" % ("H", Rx, Ry, Rz))
 
@@ -1238,7 +1238,7 @@
 
       # debug printout
       if False:
-        for key in position['template'].keys():
+        for key in list(position['template'].keys()):
           Rx, Ry, Rz = position['template'][key]
           pka_print("ATOM    669  %-3s GLU C  89    %8.3lf%8.3lf%8.3lf" % ("H", Rx, Ry, Rz))
 
@@ -1258,12 +1258,12 @@
 
       # debug printout
       if False:
-        for key in position['template'].keys():
+        for key in list(position['template'].keys()):
           Rx, Ry, Rz = position['template'][key]
           pka_print("ATOM    669  %-3s GLU B  89    %8.3lf%8.3lf%8.3lf" % ("H", Rx, Ry, Rz))
 
       # initializing iterations: copy over 'template' to 'new' & 'old'
-      for key in position['template'].keys():
+      for key in list(position['template'].keys()):
         position['new'][key] = [position['template'][key][0], position['template'][key][1], position['template'][key][2]]
         position['old'][key] = [position['template'][key][0], position['template'][key][1], position['template'][key][2]]
 
@@ -1281,7 +1281,7 @@
         if rmsd_new < rmsd_old:
           # pka_print("REMARK   update structure (iter %3d (%8.3lf%8.3lf))" % (iter, rmsd_new, rmsd_old))
           chainNumb += 1
-          for key in position['old'].keys():
+          for key in list(position['old'].keys()):
             for i in range(3):
               position['old'][key][i] = position['new'][key][i]
             if False:
@@ -1292,7 +1292,7 @@
           copied_residue.translate(dR)
         else:
           # pka_print("no update        (iter %3d (%8.3lf%8.3lf))" % (iter, rmsd_new, rmsd_old))
-          for key in position['new'].keys():
+          for key in list(position['new'].keys()):
             for i in range(3):
               position['new'][key][i] = position['old'][key][i]
 
--- a/propka30/Source/output.py
+++ b/propka30/Source/output.py
@@ -37,7 +37,7 @@
 #   Journal of Chemical Theory and Computation, 7, 525-537 (2011)
 #-------------------------------------------------------------------------------------------------------
 import sys
-import lib
+from . import lib
 revision = 182
 
 
@@ -63,7 +63,7 @@
         filename = "%s.pdb" % (protein.name)
       file = open(filename, 'w')
       if options.verbose:
-          print("writing pdbfile %s" % (filename))
+          print(("writing pdbfile %s" % (filename)))
       close_file = True
     else:
       # don't close the file, it was opened in a different place
@@ -156,7 +156,7 @@
       filename = "%s.pka" % (protein.name)
     file = open(filename, 'w')
     if verbose == True:
-      print("writing pkafile %s" % (filename))
+      print(("writing pkafile %s" % (filename)))
 
     # writing propka header
     str  = "%s\n" % ( getPropkaHeader() )
@@ -213,7 +213,7 @@
     # geting the determinants section
     str = getDeterminantSection(protein)
     if verbose:
-        print(str[:-1])
+        print((str[:-1]))
 
     str = getSummarySection(protein)
     if verbose:
--- a/propka30/Source/parameters.py
+++ b/propka30/Source/parameters.py
@@ -37,10 +37,10 @@
 #   Journal of Chemical Theory and Computation, 7, 525-537 (2011)
 #-------------------------------------------------------------------------------------------------------
 import math
-import lib
+from . import lib
 pka_print = lib.pka_print
 import sys, os
-from calculator import calculate
+from .calculator import calculate
 
 
 def pKa_mod(resName):
@@ -668,7 +668,7 @@
                    'TRP': "TRP"}
 
 
-        if resName in resType.keys():
+        if resName in list(resType.keys()):
             return resType[resName]
         else:
             return None
@@ -1019,7 +1019,7 @@
         
         if resName in Q:
             return Q[resName]
-        elif resName in self.ions.keys():
+        elif resName in list(self.ions.keys()):
             return self.ions[resName]
         else:
             return 0.00
@@ -1046,9 +1046,9 @@
                    'GLN': "AMD",
                    'TRP': "TRP"}
 
-        if resName in resType.keys():
+        if resName in list(resType.keys()):
             return resType[resName]
-        elif resName in self.ions.keys():
+        elif resName in list(self.ions.keys()):
             return 'LIG'
         else:
             return None
--- a/propka30/Source/parameters_new.py
+++ b/propka30/Source/parameters_new.py
@@ -37,7 +37,7 @@
 #   Journal of Chemical Theory and Computation, 7, 525-537 (2011)
 #-------------------------------------------------------------------------------------------------------
 
-from lib import pka_print
+from .lib import pka_print
 
 def resName2Type(resName=None):
     """
@@ -368,7 +368,7 @@
 
 
             # updating parameter matrix to full matrix
-            keys = parameters.keys()
+            keys = list(parameters.keys())
             for key1 in keys:
               for key2 in keys:
                 if key2 not in parameters[key1]:
--- a/propka30/Source/parameters_std.py
+++ b/propka30/Source/parameters_std.py
@@ -37,7 +37,7 @@
 #   Journal of Chemical Theory and Computation, 7, 525-537 (2011)
 #-------------------------------------------------------------------------------------------------------
 
-from lib import pka_print
+from .lib import pka_print
 
 def resName2Type(resName=None):
     """
@@ -362,7 +362,7 @@
 
 
             # updating side-chain parameter matrix to full matrix
-            keys = parameters.keys()
+            keys = list(parameters.keys())
             for key1 in keys:
               for key2 in keys:
                 if key2 not in parameters[key1]:
--- a/propka30/Source/pdb.py
+++ b/propka30/Source/pdb.py
@@ -37,7 +37,7 @@
 #   Journal of Chemical Theory and Computation, 7, 525-537 (2011)
 #-------------------------------------------------------------------------------------------------------
 import string, sys, copy
-import lib
+from . import lib
 pka_print = lib.pka_print
 
 excluded_resNames = ["H2O", "HOH", "SO4", "PO4", "PEG", "EPE", "NAG", "TRS"]
@@ -109,8 +109,8 @@
       for chainID in sorted( atoms.keys() ):
         for key in atoms[chainID]["keys"]:
           for atom in atoms[chainID][key]:
-            str = "%s%4d  %4s%3d%7s" % (atom.resName, atom.resNumb, atom.name, len(atom.configurations.keys()), atom.type)
-            for key in atom.configurations.keys():
+            str = "%s%4d  %4s%3d%7s" % (atom.resName, atom.resNumb, atom.name, len(list(atom.configurations.keys())), atom.type)
+            for key in list(atom.configurations.keys()):
               str += "%5s" % (key)
             pka_print(str)
 
@@ -361,7 +361,7 @@
             self.x += vector[0]
             self.y += vector[1]
             self.z += vector[2]
-            for key in self.configurations.keys():
+            for key in list(self.configurations.keys()):
               for i in range(3):
                 self.configurations[key][i] += vector[i]
 
@@ -405,14 +405,14 @@
               self.z = self.configurations[key][2]
             elif len(self.configurations) == 1:
               # get single key if only one configuration: saving back-bone protonation when previous residue doesn't have 'key'
-              for default_key in self.configurations.keys():
+              for default_key in list(self.configurations.keys()):
                 break
               self.x = self.configurations[default_key][0]
               self.y = self.configurations[default_key][1]
               self.z = self.configurations[default_key][2]
             elif True:
               # get single key if only one configuration: saving back-bone protonation when previous residue doesn't have 'key'
-              for default_key in self.configurations.keys():
+              for default_key in list(self.configurations.keys()):
                 if key[:-2] == default_key[:-2]:
                   break
               self.x = self.configurations[default_key][0]
@@ -421,7 +421,7 @@
             else:
               resLabel = "%-3s%4d%2s" % (self.resName, self.resNumb, self.chainID)
               keys = ""
-              for item in self.configurations.keys(): keys += "%5s" % (item)
+              for item in list(self.configurations.keys()): keys += "%5s" % (item)
               pka_print("configuration '%s' not found in '%s' atom '%s' [%s]" % (key, resLabel, self.name, keys))
               sys.exit(8)
 
--- a/propka30/Source/protein.py
+++ b/propka30/Source/protein.py
@@ -37,15 +37,15 @@
 #   Journal of Chemical Theory and Computation, 7, 525-537 (2011)
 #-------------------------------------------------------------------------------------------------------
 import math, sys, os, time, string
-import lib
+from . import lib
 pka_print = lib.pka_print
-import determinants
-import pdb 
+from . import determinants
+from . import pdb 
 #import debug
-import output
-import coupled_residues
-import calculator as calculate
-from chain import Chain
+from . import output
+from . import coupled_residues
+from . import calculator as calculate
+from .chain import Chain
 
 
 
@@ -290,7 +290,7 @@
         """ 
         protonates the protein according to given scheme
         """ 
-        from protonator import makeProtonator
+        from .protonator import makeProtonator
         please = makeProtonator(scheme=scheme)
         self.removeHydrogens()
 
@@ -317,7 +317,7 @@
         """ 
         # create a default version if not provided
         if version == None:
-          import version
+          from . import version
           version = version.makeVersion(label=options.version_label)
 
         if len(self.configurations) == 1:
@@ -408,25 +408,25 @@
           pka_print(str)
 
         # get the average pKa properties by dividing the sum with number of configurations, len(configuration keys)
-        from determinants import Determinant
+        from .determinants import Determinant
         for residue in pka_residues:
           residue_key = residue.label
           sum_pka = 0.00; sum_Nmass = 0.00; sum_Emass = 0.00; sum_Nlocl = 0.00; sum_Elocl = 0.00
-          for key in pkas[residue_key].keys():
+          for key in list(pkas[residue_key].keys()):
             sum_pka   += pkas[residue_key][key]
             sum_Nmass += Nmass[residue_key][key]
             sum_Emass += Emass[residue_key][key]
             sum_Nlocl += Nlocl[residue_key][key]
             sum_Elocl += Elocl[residue_key][key]
-          residue.pKa_pro = sum_pka/len( pkas[residue_key].keys() )
-          residue.Nmass   = sum_Nmass/len( pkas[residue_key].keys() )
-          residue.Emass   = sum_Emass/len( pkas[residue_key].keys() )
-          residue.Nlocl   = sum_Nlocl/len( pkas[residue_key].keys() )
-          residue.Elocl   = sum_Elocl/len( pkas[residue_key].keys() )
+          residue.pKa_pro = sum_pka/len( list(pkas[residue_key].keys()) )
+          residue.Nmass   = sum_Nmass/len( list(pkas[residue_key].keys()) )
+          residue.Emass   = sum_Emass/len( list(pkas[residue_key].keys()) )
+          residue.Nlocl   = sum_Nlocl/len( list(pkas[residue_key].keys()) )
+          residue.Elocl   = sum_Elocl/len( list(pkas[residue_key].keys()) )
           residue.determinants = [[], [], []]
           for type in range(3):
-            for key in determinants[residue_key][type].keys():
-              value = determinants[residue_key][type][key] / len( pkas[residue_key].keys() )
+            for key in list(determinants[residue_key][type].keys()):
+              value = determinants[residue_key][type][key] / len( list(pkas[residue_key].keys()) )
               if abs(value) > 0.005:  # <-- removing determinant that appears as 0.00
                 newDeterminant = Determinant(key, value)
                 residue.determinants[type].append(newDeterminant)
@@ -774,7 +774,7 @@
         """ 
         mutates the protein according to 'mutation' using 'method'
         """
-        import mutate
+        from . import mutate
         newProtein = mutate.makeMutatedProtein(self, mutation=mutation, atoms=atoms, options=options)
 
         return  newProtein
@@ -784,7 +784,7 @@
         """ 
         permutes multiple mutations and determins the most stable combination; note, you need the version for stability calculations
         """
-        import mutate
+        from . import mutate
         best_mutation = mutate.optimizeMutationDeterminants(self, mutation=mutation, atoms=atoms, alignment=alignment, version=version, options=options)
 
         return  best_mutation
@@ -794,7 +794,7 @@
         """ 
         permutes multiple mutations and determins the most stable combination; note, you need the version for stability calculations
         """
-        import mutate
+        from . import mutate
         best_mutation = mutate.optimizeMultipleMutations(self, mutations=mutations, atoms=atoms, alignment=None, version=version, options=options)
 
         return  best_mutation
@@ -914,7 +914,7 @@
         file.write(str)
 
         # set 'protein' based on pdbcode
-        for protein in experiment.keys():
+        for protein in list(experiment.keys()):
           str  = "%s" % (protein)
           #pka_print(str)
           #str += " %s" % (experiment[protein]['pdb'])
@@ -925,7 +925,7 @@
         # setting a list of labels to work with
         if   labels == None or labels == "ALL":
           labels = []
-          for label in experiment[protein].keys():
+          for label in list(experiment[protein].keys()):
             labels.append(label)
 
         # iterating that list
@@ -961,7 +961,7 @@
     """ 
     Reads necessary information about residues (includes ions)
     """ 
-    from parameters_new import resName2Type, getQs, pKa_mod
+    from .parameters_new import resName2Type, getQs, pKa_mod
     resInfo = {}
     # reading residue information from parameters.py
     resInfo['resType'] = resName2Type()
@@ -970,7 +970,7 @@
     resInfo['type']    = {'C- ': "C-terminus",
                           'N+ ': "N-terminus"}
     # setting up 'type' = 'amino-acid' for residues
-    for resName in resInfo['resType'].keys():
+    for resName in list(resInfo['resType'].keys()):
       if resName not in resInfo['type']:
         resInfo['type'][resName] = "amino-acid"
 
--- a/propka30/Source/protonate.py
+++ b/propka30/Source/protonate.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 #
 # * This library is free software; you can redistribute it and/or
 # * modify it under the terms of the GNU Lesser General Public
@@ -39,9 +39,9 @@
 #-------------------------------------------------------------------------------------------------------
 
 
-from vector_algebra import *
-import bonds, pdb
-from lib import pka_print
+from .vector_algebra import *
+from . import bonds, pdb
+from .lib import pka_print
 
 class Protonate:
     """ Protonates atoms using VSEPR theory """
--- a/propka30/Source/protonator.py
+++ b/propka30/Source/protonator.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 #
 # * This library is free software; you can redistribute it and/or
 # * modify it under the terms of the GNU Lesser General Public
@@ -40,11 +40,11 @@
 
 
 import  sys, math
-from vector_algebra import *
-import bonds as bonds
-import pdb as pdb
+from .vector_algebra import *
+from . import bonds as bonds
+from . import pdb as pdb
 
-from lib import pka_print
+from .lib import pka_print
 
 
 def makeProtonator(scheme=None):
@@ -201,7 +201,7 @@
         X1, X2, X3 = atoms
         H  = X1.makeCopy(name=name, element='H')
 
-        for key in H.configurations.keys():
+        for key in list(H.configurations.keys()):
           for atom in [H, X1, X2, X3]:
             atom.setConfiguration(key)
 
@@ -227,7 +227,7 @@
         X1, X2, X3 = atoms
         H  = X1.makeCopy(name=name, element='H')
 
-        for key in H.configurations.keys():
+        for key in list(H.configurations.keys()):
           for atom in [H, X1, X2, X3]:
             atom.setConfiguration(key)
 
@@ -253,7 +253,7 @@
         X1, X2, X3 = atoms
         H  = X2.makeCopy(name=name, element='H')
 
-        for key in H.configurations.keys():
+        for key in list(H.configurations.keys()):
           for atom in [H, X1, X2, X3]:
             atom.setConfiguration(key)
 
--- a/propka30/Source/residue.py
+++ b/propka30/Source/residue.py
@@ -40,8 +40,8 @@
 import string
 import math
 import copy
-import lib
-from pdb import Atom
+from . import lib
+from .pdb import Atom
 pka_print = lib.pka_print
 
 class Residue:
@@ -107,7 +107,7 @@
                 self.type = "ligand"
 
             # setting the number of configurations for this residue
-            for key in atom.configurations.keys():
+            for key in list(atom.configurations.keys()):
                 if key not in self.configurations:
                     self.configurations.append(key)
 
@@ -792,7 +792,7 @@
           atom.x += translation[0]
           atom.y += translation[1]
           atom.z += translation[2]
-          for key in atom.configurations.keys():
+          for key in list(atom.configurations.keys()):
             for i in range(3):
               atom.configurations[key][i] += translation[i]
 
@@ -801,7 +801,7 @@
         """
         rotate residue theta radians around axis with center=center
         """
-        from rotate import generalRotationMatrix
+        from .rotate import generalRotationMatrix
         translate = [0.00, 0.00, 0.00]
         number_of_atoms = 0
         for atom in self.atoms:
@@ -839,7 +839,7 @@
           atom.z = new_position[2]
 
           # rotate configuration
-          for key in atom.configurations.keys():
+          for key in list(atom.configurations.keys()):
             for xyz in range(3):
               new_position[xyz] = translate[xyz]
               for i in range(3):
@@ -855,7 +855,7 @@
         """
         making a copy of this residue
         """
-        from protein import getResidueParameters
+        from .protein import getResidueParameters
         if chainID == None: chainID = self.chainID
         if resNumb == None: resNumb = self.resNumb
 
--- a/propka30/Source/rotate.py
+++ b/propka30/Source/rotate.py
@@ -40,7 +40,7 @@
 
 
 import math, os, sys, random
-from lib import pka_print
+from .lib import pka_print
 
 
 def generateCorrespondingAtoms():
@@ -183,7 +183,7 @@
     """
     translates the position according to 'translation'
     """
-    for key in position.keys():
+    for key in list(position.keys()):
       for i in range(3):
         position[key][i] += translation[i]
 
@@ -200,7 +200,7 @@
         translate[i] += position[key][i]/len(center)
 
     # translate to rotation center
-    for key in position.keys():
+    for key in list(position.keys()):
       for i in range(3):
         position[key][i] -= translate[i]
 
@@ -209,7 +209,7 @@
 
     # do the actual rotation
     new_position = [None, None, None]
-    for key in position.keys():
+    for key in list(position.keys()):
       # rotate
       for xyz in range(3):
         new_position[xyz] = translate[xyz]
@@ -270,7 +270,7 @@
       atom.z = new_position[2]
 
       # rotate configuration
-      for key in atom.configurations.keys():
+      for key in list(atom.configurations.keys()):
         for xyz in range(3):
           new_position[xyz] = translate[xyz]
           for i in range(3):
--- a/propka30/Source/vector_algebra.py
+++ b/propka30/Source/vector_algebra.py
@@ -37,7 +37,7 @@
 #   Journal of Chemical Theory and Computation, 7, 525-537 (2011)
 #-------------------------------------------------------------------------------------------------------
 import math, sys
-import lib
+from . import lib
 pka_print = lib.pka_print
 
 
@@ -275,9 +275,9 @@
 
         # store vectors for all configurations of atoms
         if atom1!=0:
-            self.keys = lib.get_sorted_configurations(atom1.configurations.keys())
+            self.keys = lib.get_sorted_configurations(list(atom1.configurations.keys()))
             if atom2!=0:
-                keys2 = lib.get_sorted_configurations(atom2.configurations.keys())
+                keys2 = lib.get_sorted_configurations(list(atom2.configurations.keys()))
                 if self.keys != keys2:
                     pka_print("ERROR: Inequivalent configurations for atoms, please correct your pdbfile to single configuration")
                     pka_print("%s\n%s" % (atom1, atom2))
--- a/propka30/Source/version.py
+++ b/propka30/Source/version.py
@@ -38,9 +38,9 @@
 #-------------------------------------------------------------------------------------------------------
 
 import math
-import lib
+from . import lib
 import sys, os
-import calculator as calculate
+from . import calculator as calculate
 pka_print = lib.pka_print
 
 
@@ -413,7 +413,7 @@
         if resType == None:
           # doing it for all residues
           #pka_print("changing back-bone parameters")
-          for key in self.BackBoneParameters.keys():
+          for key in list(self.BackBoneParameters.keys()):
             self.BackBoneParameters[key][0] = dpka
             str  = " %s: %6.2lf " % (key, self.BackBoneParameters[key][0])
             str += "[%5.2lf,%5.2lf]" % (self.BackBoneParameters[key][1][0], self.BackBoneParameters[key][1][1])
@@ -434,9 +434,9 @@
         """
         if resType == None:
           # doing it for all residues
-          for key1 in self.SideChainParameters.keys():
+          for key1 in list(self.SideChainParameters.keys()):
             #pka_print("changing side-chain parameters for resType \"%s\"" % (key1))
-            for key2 in self.SideChainParameters[key1].keys():
+            for key2 in list(self.SideChainParameters[key1].keys()):
               self.SideChainParameters[key1][key2][0] = dpka
               str  = " %s: %6.2lf " % (key2, self.SideChainParameters[key1][key2][0])
               str += "[%5.2lf,%5.2lf]" % (self.SideChainParameters[key1][key2][1][0], self.SideChainParameters[key1][key2][1][1])
@@ -445,7 +445,7 @@
           # doing it just for residue 'key1'
           key1 = resType
           #pka_print("changing side-chain parameters for resType \"%s\"" % (key1))
-          for key2 in self.SideChainParameters[key1].keys():
+          for key2 in list(self.SideChainParameters[key1].keys()):
             self.SideChainParameters[key1][key2][0] = dpka
             self.SideChainParameters[key2][key1][0] = dpka
             str  = " %s: %6.2lf " % (key2, self.SideChainParameters[key1][key2][0])
@@ -457,7 +457,7 @@
         """
         setting the Coulomb model, and its parameters. If parameter is not defined, it will take the default from 'parameters.py'
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
         coulomb_parameters = parameters.getCoulombParameters()
         if label not in coulomb_parameters:
           pka_print("do not accept Coulomb model \"%s\\n" % (label))
@@ -494,7 +494,7 @@
         """
         setting the desolvation model, and its parameters
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
         desolvation_parameters = parameters.getDesolvationParameters()
         if label not in desolvation_parameters:
           pka_print("do not accept solvation model \"%s\\n" % (label))
@@ -661,7 +661,7 @@
         """
         Rules of action for version Jan01
         """
-        import parameters_std as parameters
+        from . import parameters_std as parameters
     
         self.name             = "Jan01"
         self.Nmin             =  300
@@ -750,7 +750,7 @@
         """
         Rules of action for version Jan15
         """
-        import parameters_std as parameters
+        from . import parameters_std as parameters
 
         self.name             = "Jan15"
         self.Nmin             =  300
@@ -839,7 +839,7 @@
         """
         Rules of action for version May13
         """
-        import parameters_std as parameters
+        from . import parameters_std as parameters
 
         self.name             = "May13"
         self.coulomb_cutoff =   7.00
@@ -869,7 +869,7 @@
         """
         Rules of action for version Dec18
         """
-        import parameters_std as parameters
+        from . import parameters_std as parameters
 
         self.name             = "Dec18"
         self.Nmin             =  300
@@ -925,7 +925,7 @@
         """
         Rules of action for version Dec19
         """
-        import parameters_std as parameters
+        from . import parameters_std as parameters
         
         self.name             = "Dec19"
         self.Nmin             =  300
@@ -981,7 +981,7 @@
         """
         Rules of action for version Aug24
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Aug24"
         self.Nmin             =  280
@@ -1009,7 +1009,7 @@
         """
         Rules of action for version Aug30
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Aug30"
         self.Nmin             =  280
@@ -1037,7 +1037,7 @@
         """
         Rules of action for version Aug31
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Aug31"
         self.Nmin             =  280
@@ -1065,7 +1065,7 @@
         """
         Rules of action for version Sep05
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Sep05"
         self.Nmin             =  280
@@ -1093,7 +1093,7 @@
         """
         Rules of action for version Sep06
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Sep06"
         self.Nmin             =  280
@@ -1121,7 +1121,7 @@
         """
         Rules of action for version Sep07
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Sep07"
         self.Nmin             =  280
@@ -1149,7 +1149,7 @@
         """
         Rules of action for version Sep08
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Sep08"
         self.Nmin             =  280
@@ -1177,7 +1177,7 @@
         """
         Rules of action for version Oct13, based on Sep07
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Oct13"
         self.Nmin             =  280
@@ -1239,7 +1239,7 @@
         
         if resName in Q:
             return Q[resName]
-        elif resName in self.ions.keys():
+        elif resName in list(self.ions.keys()):
             return self.ions[resName]
         else:
             return 0.00
@@ -1266,9 +1266,9 @@
                    'GLN': "AMD",
                    'TRP': "TRP"}
 
-        if resName in resType.keys():
+        if resName in list(resType.keys()):
             return resType[resName]
-        elif resName in self.ions.keys():
+        elif resName in list(self.ions.keys()):
             return 'LIG'
         else:
             return None
@@ -1309,7 +1309,7 @@
         """
         Rules of action for version Oct14
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Oct14"
         self.Nmin             =  280
@@ -1337,7 +1337,7 @@
         """
         Rules of action for version Oct14
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Nov28"
         self.Nmin             =  280
@@ -1365,7 +1365,7 @@
         """
         Rules of action for version Nov29
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Nov29"
         self.Nmin             =  280
@@ -1393,7 +1393,7 @@
         """
         Rules of action for version Nov30
         """
-        import parameters_new as parameters
+        from . import parameters_new as parameters
 
         self.name             = "Nov30"
         self.Nmin             =  280
--- a/propka30/experiment.py
+++ b/propka30/experiment.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
 #
 # * This library is free software; you can redistribute it and/or
 # * modify it under the terms of the GNU Lesser General Public
--- a/propka30/propka.py
+++ b/propka30/propka.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 #
 # * This library is free software; you can redistribute it and/or
 # * modify it under the terms of the GNU Lesser General Public
@@ -41,7 +41,7 @@
 
 import string, re, sys, os, math
 
-import Source.lib as lib
+from Source import lib
 from Source.protein import Protein
  
 
--- a/propka30/propka_det.py
+++ b/propka30/propka_det.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
 #
 # * This library is free software; you can redistribute it and/or
 # * modify it under the terms of the GNU Lesser General Public
@@ -39,8 +39,8 @@
 #-------------------------------------------------------------------------------------------------------
 
 import string, re, sys, os, math
-import Source.version as propka
-import Source.lib as lib
+from Source import version as propka
+from Source import lib
 from Source.protein import Protein
 from Source.mutate import makeCompositeAtomsDictionary
  
--- a/propka30/propka_gui.py
+++ b/propka30/propka_gui.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
 #
 # * This library is free software; you can redistribute it and/or
 # * modify it under the terms of the GNU Lesser General Public
@@ -39,8 +39,8 @@
 #-------------------------------------------------------------------------------------------------------
 
 import string, re, sys, os, math
-import Source.version as propka
-import Source.lib as lib
+from Source import version as propka
+from Source import lib
 from Source.protein import Protein
 from Source.mutate import makeCompositeAtomsDictionary
 pka_print = lib.pka_print
--- a/propka30/propka_mut.py
+++ b/propka30/propka_mut.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
 #
 # * This library is free software; you can redistribute it and/or
 # * modify it under the terms of the GNU Lesser General Public
@@ -40,11 +40,11 @@
 
 import string, re, sys, os, math
 
-import Source.lib as lib
-import Source.mutate as mutate
+from Source import lib
+from Source import mutate
 from Source.protein import Protein
 from Source.pdb import readPDB
-import Source.version as propka
+from Source import version as propka
 pka_print = lib.pka_print
 
 
--- a/propka30/propka_pdb.py
+++ b/propka30/propka_pdb.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/python3
 #
 # * This library is free software; you can redistribute it and/or
 # * modify it under the terms of the GNU Lesser General Public
@@ -40,7 +40,7 @@
 
 import string, re, sys, os, math
 
-import Source.lib as lib
+from Source import lib
 from Source.protein import Protein
  
 
--- a/querystatus.py
+++ b/querystatus.py
@@ -9,7 +9,7 @@
 import sys
 import cgi
 import cgitb
-import os,shutil,glob,string,time,urllib
+import os,shutil,glob,string,time,urllib.request,urllib.parse,urllib.error
 from datetime import timedelta
 from src.server import *
 from src.aconf import *
@@ -72,7 +72,7 @@
         out_f.write(origin_line)
 
 
-        for x in xrange(3):
+        for x in range(3):
             split_line = in_f.readline().split()
             grid_dims = [float(item) for item in split_line[-3:]]
 
@@ -146,7 +146,7 @@
         # construct soap request
         try:
             status=appServicePort.queryStatus(queryStatusRequest(jobid))
-        except Exception, e:
+        except Exception as e:
             return ["error"]
         if status._code == 4:
             return ["error"]
@@ -170,16 +170,16 @@
         Main method for determining the query page output
     """
     logopts = {}
-    print "Content-type: text/html\n\n"
+    print("Content-type: text/html\n\n")
     calctype = form["calctype"].value
 
     # prints version error, if it exists
     if form["jobid"].value == 'False':
-        print printheader("%s Job Status Page" % calctype.upper())
+        print(printheader("%s Job Status Page" % calctype.upper()))
         progress = "version_mismatch"
         runtime = 0
     elif form["jobid"].value == 'notenoughmem':
-        print printheader("%s Job Status Page" % calctype.upper())
+        print(printheader("%s Job Status Page" % calctype.upper()))
         progress = "not_enough_memory"
         runtime = 0
     else:
@@ -194,7 +194,7 @@
         string+= "\t\t<meta http-equiv=\"Refresh\" content=\"0; url=%s%s%s.html\">\n" % (WEBSITE, TMPDIR, form["jobid"].value)
         string+= "\t</head>\n"
         string+= "</html>\n"
-        print string
+        print(string)
         return
 
     # prepares for Opal query, if necessary
@@ -283,18 +283,18 @@
             resultsurl = '%squerystatus.cgi?jobid=%s&calctype=apbs' % (WEBSITE, form["jobid"].value)
 
     if progress == "complete":
-        print printheader("%s Job Status Page" % calctype.upper(), jobid=form["jobid"].value)
+        print(printheader("%s Job Status Page" % calctype.upper(), jobid=form["jobid"].value))
 
     elif progress == "error":
-        print printheader("%s Job Status Page - Error" % calctype.upper(),0, jobid=form["jobid"].value)
+        print(printheader("%s Job Status Page - Error" % calctype.upper(),0, jobid=form["jobid"].value))
 
     elif progress == "running": # job is not complete, refresh in 30 seconds
-        print printheader("%s Job Status Page" % calctype.upper(), refresh, jobid=form["jobid"].value)
+        print(printheader("%s Job Status Page" % calctype.upper(), refresh, jobid=form["jobid"].value))
 
-    print "<BODY>\n<P>"
-    print "<p></p>"
-    print '<div id="content">'
-    print "<h3>Status:"
+    print("<BODY>\n<P>")
+    print("<p></p>")
+    print('<div id="content">')
+    print("<h3>Status:")
 
     color = "FA3434"
     image = WEBSITE+"images/red_x.png"
@@ -306,12 +306,12 @@
         color = "ffcc00"
         image = WEBSITE+"images/yellow_exclamation.png"
 
-    print "<strong style=\"color:#%s;\">%s</strong>" % (color, progress)
-    print "<img src=\"%s\"><br />" % image
-    print "</h3>"
-    print "Run time: " + str(timedelta(seconds=round(runtime))) + '<br />'
-    print "Current time: %s<br />" % time.asctime()
-    print "</P>\n<HR>\n<P>"
+    print("<strong style=\"color:#%s;\">%s</strong>" % (color, progress))
+    print("<img src=\"%s\"><br />" % image)
+    print("</h3>")
+    print("Run time: " + str(timedelta(seconds=round(runtime))) + '<br />')
+    print("Current time: %s<br />" % time.asctime())
+    print("</P>\n<HR>\n<P>")
 
     if progress == "complete":
         if calctype=="pdb2pqr":
@@ -325,8 +325,8 @@
             resp = appServicePort.getOutputs(getOutputsRequest(jobid))
             filelist = resp._outputFile
 
-        print "Here are the results:<ul>"
-        print "<li>Input files<ul>"
+        print("Here are the results:<ul>")
+        print("<li>Input files<ul>")
 
         if calctype=="pdb2pqr":
             # this code should be cleaned up once local PDB2PQR runs output the PDB file with the .pdb extension
@@ -335,26 +335,26 @@
                     file_name = filelist[i]._name
                     if ((len(file_name) == 4 and '.' not in file_name) or
                         (file_name.lower().endswith(".pdb") and "pdb2pka_output" not in file_name)):
-                        print "<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name)
+                        print("<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name))
 
                     if file_name.lower().endswith((".mol", ".mol2", ".names", ".dat")) and "pdb2pka_output" not in file_name:
-                        print "<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name)
+                        print("<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name))
 
 
             else:
-                print "<li><a href=%s%s%s/%s.pdb>%s.pdb</a></li>" % (WEBSITE, TMPDIR, jobid, jobid, jobid)
+                print("<li><a href=%s%s%s/%s.pdb>%s.pdb</a></li>" % (WEBSITE, TMPDIR, jobid, jobid, jobid))
 
         elif calctype=="apbs":
             if have_opal:
                 for i in range(0,len(filelist)):
                     if filelist[i]._name == "apbsinput.in" or filelist[i]._name[-4:] == ".pqr":
-                        print "<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name)
+                        print("<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name))
             else:
-                print "<li><a href=%s%s%s/apbsinput.in>apbsinput.in</a></li>" % (WEBSITE, TMPDIR, jobid)
-                print "<li><a href=%s%s%s/%s.pqr>%s.pqr</a></li>" % (WEBSITE, TMPDIR, jobid, jobid, jobid)
+                print("<li><a href=%s%s%s/apbsinput.in>apbsinput.in</a></li>" % (WEBSITE, TMPDIR, jobid))
+                print("<li><a href=%s%s%s/%s.pqr>%s.pqr</a></li>" % (WEBSITE, TMPDIR, jobid, jobid, jobid))
 
-        print "</ul></li>"
-        print "<li>Output files<ul>"
+        print("</ul></li>")
+        print("<li>Output files<ul>")
 
         queryString = [str(os.environ["REMOTE_ADDR"])]
         # Getting PDB2PQR run log info
@@ -369,17 +369,17 @@
             if have_opal:
                 for i in range(0,len(filelist)):
                     if filelist[i]._name.endswith((".propka", "-typemap.html")):
-                        print "<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name)
+                        print("<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name))
 
                     if filelist[i]._name.endswith(".in") and "pdb2pka_output" not in filelist[i]._name:
-                        print "<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name)
+                        print("<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name))
 
                     if filelist[i]._name.endswith(".pqr") and not filelist[i]._name.endswith("background_input.pqr"):
-                        print "<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name)
+                        print("<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name))
 
                     #Get the first line of the summary file.
                     if filelist[i]._name.endswith(".summary"):
-                        f=urllib.urlopen(filelist[i]._url)
+                        f=urllib.request.urlopen(filelist[i]._url)
                         summaryLine = f.readline().strip()
                         #logopts["pdb"]=logopts.get("pdb", "") + "|" + summaryLine
                         queryString.append(summaryLine)
@@ -403,7 +403,7 @@
                             continue
                         outputfilelist.append('%s%s' % (jobid, extension))
                 for outputfile in outputfilelist:
-                    print "<li><a href=%s%s%s/%s>%s</a></li>" % (WEBSITE, TMPDIR, jobid, outputfile, outputfile)
+                    print("<li><a href=%s%s%s/%s>%s</a></li>" % (WEBSITE, TMPDIR, jobid, outputfile, outputfile))
 
             logopts['queryPDB2PQR'] = '|'.join(queryString)
 
@@ -417,7 +417,7 @@
                         currentpath = os.getcwd()
                         zipjobid = filelist[i]._name.split("-")[0]
                         dxfilename = '%s%s%s/%s' % (INSTALLDIR, TMPDIR, zipjobid, filelist[i]._name)
-                        urllib.urlretrieve(filelist[i]._url, dxfilename)
+                        urllib.request.urlretrieve(filelist[i]._url, dxfilename)
                         os.chdir('%s%s%s' % (INSTALLDIR, TMPDIR, zipjobid))
                         # making both the dx file and the compressed file (.gz) available in the directory
                         syscommand = 'cp %s dxbkupfile' % (filelist[i]._name)
@@ -444,8 +444,8 @@
                         os.chdir(currentpath)
                         outputcubefilezip = cubefilebasename+".gz"
 
-                        print "<li><a href=%s%s%s/%s>%s</a></li>" % (WEBSITE, TMPDIR, zipjobid, outputfilezip, outputfilezip)
-                        print "<li><a href=%s%s%s/%s>%s</a></li>" % (WEBSITE, TMPDIR, zipjobid, outputcubefilezip, outputcubefilezip)
+                        print("<li><a href=%s%s%s/%s>%s</a></li>" % (WEBSITE, TMPDIR, zipjobid, outputfilezip, outputfilezip))
+                        print("<li><a href=%s%s%s/%s>%s</a></li>" % (WEBSITE, TMPDIR, zipjobid, outputcubefilezip, outputcubefilezip))
 
             else:
                 outputfilelist = glob.glob('%s%s%s/%s-*.dx' % (INSTALLDIR, TMPDIR, jobid, jobid))
@@ -472,7 +472,7 @@
 
                     createcube(dxfile, pqrfilename, cubefilename)
 
-                    print "<li><a href=%s%s%s/%s>%s</a></li>" % (WEBSITE, TMPDIR, jobid, os.path.basename(outputfilezip), os.path.basename(outputfilezip))
+                    print("<li><a href=%s%s%s/%s>%s</a></li>" % (WEBSITE, TMPDIR, jobid, os.path.basename(outputfilezip), os.path.basename(outputfilezip)))
 
                 outputcubefilelist = glob.glob('%s%s%s/%s.cube' % (INSTALLDIR, TMPDIR, jobid, jobid))
                 for cubefile in outputcubefilelist:
@@ -489,7 +489,7 @@
                     os.chdir(currentpath)
                     outputcubefilezip = cubefile+".gz"
 
-                    print "<li><a href=%s%s%s/%s>%s</a></li>" % (WEBSITE, TMPDIR, jobid, os.path.basename(outputcubefilezip), os.path.basename(outputcubefilezip))
+                    print("<li><a href=%s%s%s/%s>%s</a></li>" % (WEBSITE, TMPDIR, jobid, os.path.basename(outputcubefilezip), os.path.basename(outputcubefilezip)))
 
             logopts['queryAPBS'] = '|'.join(queryString)
 
@@ -501,22 +501,22 @@
                         outputfilelist.append((filelist[i]._url, os.path.basename(filelist[i]._name)))
                         #print "<li><a href=%s>%s</a></li>" % (filelist[i]._url, filelist[i]._name)
                 if outputfilelist:
-                    print "</ul></li>"
-                    print "<li>PDB2PKA files<ul>"
+                    print("</ul></li>")
+                    print("<li>PDB2PKA files<ul>")
                     for outputfile in outputfilelist:
-                        print "<li><a href=%s>%s</a></li>" % (outputfile[0], outputfile[1])
+                        print("<li><a href=%s>%s</a></li>" % (outputfile[0], outputfile[1]))
             else:
                 outputfilelist = glob.glob('%s%s%s/pdb2pka_output/*.DAT' % (INSTALLDIR, TMPDIR, jobid))
                 outputfilelist.extend(glob.glob('%s%s%s/pdb2pka_output/*.txt' % (INSTALLDIR, TMPDIR, jobid)))
                 outputfilelist = [os.path.basename(outputfile) for outputfile in outputfilelist]
                 if outputfilelist:
-                    print "</ul></li>"
-                    print "<li>PDB2PKA files<ul>"
+                    print("</ul></li>")
+                    print("<li>PDB2PKA files<ul>")
                     for outputfile in outputfilelist:
-                        print "<li><a href=%s%s%s/pdb2pka_output/%s>%s</a></li>" % (WEBSITE, TMPDIR, jobid, outputfile, outputfile)
+                        print("<li><a href=%s%s%s/pdb2pka_output/%s>%s</a></li>" % (WEBSITE, TMPDIR, jobid, outputfile, outputfile))
 
-        print "</ul></li>"
-        print "<li>Runtime and debugging information<ul>"
+        print("</ul></li>")
+        print("<li>Runtime and debugging information<ul>")
 
         if have_opal:
             stdouturl = resp._stdOut
@@ -525,11 +525,11 @@
             stdouturl = "%s%s%s/%s_stdout.txt" % (WEBSITE, TMPDIR, jobid, calctype)
             stderrurl = "%s%s%s/%s_stderr.txt" % (WEBSITE, TMPDIR, jobid, calctype)
 
-        print "<li><a href=%s>Program output (stdout)</a></li>" % stdouturl
-        print "<li><a href=%s>Program errors and warnings (stderr)</a></li>" % stderrurl
+        print("<li><a href=%s>Program output (stdout)</a></li>" % stdouturl)
+        print("<li><a href=%s>Program errors and warnings (stderr)</a></li>" % stderrurl)
 
 
-        print "</ul></li></ul>"
+        print("</ul></li></ul>")
 
 
         #if have_opal:
@@ -552,17 +552,17 @@
         #            print "<li><a href=%s>%s</a></li>" % (WEBSITE+TMPDIR+jobid+"/"+line,printname)
 
         if calctype=="pdb2pqr" and apbs_input and HAVE_APBS:
-            print "</ul></p><hr><p><b><a href=%s>Click here</a> to run APBS with your results.</b></p>" % nexturl
+            print("</ul></p><hr><p><b><a href=%s>Click here</a> to run APBS with your results.</b></p>" % nexturl)
         elif calctype=="apbs":
             #print "</ul></p><hr><p><b><a href=%s>Click here</a> to visualize your results in Jmol.</b></p>" % nexturl
-            print "</ul></p><hr><p><b>Visualize your results online:"
-            print "<ul> <li><a href=%s>3Dmol</a></li><li><a href=%s>Jmol</a></li></ul>" % (url_3dmol, url_jmol)
+            print("</ul></p><hr><p><b>Visualize your results online:")
+            print("<ul> <li><a href=%s>3Dmol</a></li><li><a href=%s>Jmol</a></li></ul>" % (url_3dmol, url_jmol))
 
     elif progress == "error":
-        print "There was an error with your query request. This page will not refresh."
+        print("There was an error with your query request. This page will not refresh.")
 
-        print "</ul></li>"
-        print "<li>Runtime and debugging information<ul>"
+        print("</ul></li>")
+        print("<li>Runtime and debugging information<ul>")
 
         if have_opal:
             resp = appServicePort.getOutputs(getOutputsRequest(jobid))
@@ -573,52 +573,52 @@
             stdouturl = "%s%s%s/%s_stdout.txt" % (WEBSITE, TMPDIR, jobid, calctype)
             stderrurl = "%s%s%s/%s_stderr.txt" % (WEBSITE, TMPDIR, jobid, calctype)
 
-        print "<li><a href=%s>Program output (stdout)</a></li>" % stdouturl
-        print "<li><a href=%s>Program errors and warnings (stderr)</a></li>" % stderrurl
+        print("<li><a href=%s>Program output (stdout)</a></li>" % stdouturl)
+        print("<li><a href=%s>Program errors and warnings (stderr)</a></li>" % stderrurl)
 
-        print "</ul></li></ul>"
+        print("</ul></li></ul>")
 
         if have_opal:
-            print " <br />If your job has been running for a prolonged period of time and failed with no reason listed in the standard out or standard error, then the job probably timed out and was terminated by the system.<br />"
+            print(" <br />If your job has been running for a prolonged period of time and failed with no reason listed in the standard out or standard error, then the job probably timed out and was terminated by the system.<br />")
 
-        print '<br />If you are having trouble running PDB2PQR on the webserver, please download the <a href="http://www.poissonboltzmann.org/docs/downloads/">command line version of PDB2PQR</a> and run the job from there.'
+        print('<br />If you are having trouble running PDB2PQR on the webserver, please download the <a href="http://www.poissonboltzmann.org/docs/downloads/">command line version of PDB2PQR</a> and run the job from there.')
 
 
 
     elif progress == "running":
-        print "Page will refresh in %d seconds<br />" % refresh
-        print "<HR>"
+        print("Page will refresh in %d seconds<br />" % refresh)
+        print("<HR>")
 
         if not have_opal:
-            print "</ul></li>"
-            print "<li>Runtime and debugging information<ul>"
+            print("</ul></li>")
+            print("<li>Runtime and debugging information<ul>")
             stdouturl = "%s%s%s/%s_stdout.txt" % (WEBSITE, TMPDIR, jobid, calctype)
             stderrurl = "%s%s%s/%s_stderr.txt" % (WEBSITE, TMPDIR, jobid, calctype)
-            print "<li><a href=%s>Program output (stdout)</a></li>" % stdouturl
-            print "<li><a href=%s>Program errors and warnings (stderr)</a></li>" % stderrurl
-            print "</ul></li></ul>"
+            print("<li><a href=%s>Program output (stdout)</a></li>" % stdouturl)
+            print("<li><a href=%s>Program errors and warnings (stderr)</a></li>" % stderrurl)
+            print("</ul></li></ul>")
 
-        print "<small>Your results will appear at <a href=%s>this page</a>. If you want, you can bookmark it and come back later (note: results are only stored for approximately 12-24 hours).</small>" % resultsurl
+        print("<small>Your results will appear at <a href=%s>this page</a>. If you want, you can bookmark it and come back later (note: results are only stored for approximately 12-24 hours).</small>" % resultsurl)
     elif progress == "version_mismatch":
-        print "The versions of APBS on the local server and on the Opal server do not match, so the calculation could not be completed"
+        print("The versions of APBS on the local server and on the Opal server do not match, so the calculation could not be completed")
 
-    print "</P>"
-    print "<script type=\"text/javascript\">"
+    print("</P>")
+    print("<script type=\"text/javascript\">")
     for key in logopts:
-        print getEventTrackingString('queryData', key, logopts[key]),
+        print(getEventTrackingString('queryData', key, logopts[key]), end=' ')
         #print "_gaq.push(['_trackPageview', '/main_cgi/has_%s_%s.html']);" % (key, logopts[key])
-    print "</script>"
-    print "</div> <!--end content div-->"
-    print "</BODY>"
-    print "</HTML>"
+    print("</script>")
+    print("</div> <!--end content div-->")
+    print("</BODY>")
+    print("</HTML>")
 
-if __name__ == "__main__" and os.environ.has_key("REQUEST_METHOD"):
+if __name__ == "__main__" and "REQUEST_METHOD" in os.environ:
     """ Determine if called from command line or CGI """
     refresh=30
 
-    if not form.has_key("jobid") and form["calctype"].value=="pdb2pqr":
-        print "Content-type: text/html\n\n"
-        print printheader("PDB2PQR Job Status - Error")
+    if "jobid" not in form and form["calctype"].value=="pdb2pqr":
+        print("Content-type: text/html\n\n")
+        print(printheader("PDB2PQR Job Status - Error"))
         text="<BODY>\n"
         text+="\t<H2>Missing jobid field</H2>\n"
         text+="\t<P>Your request url is missing the jobid field</P>\n"
@@ -632,7 +632,7 @@
         text += "pageTracker._trackPageview();"
         text += "} catch(err) {}</script>"
         text+="</BODY>\n</HTML>"
-        print text
+        print(text)
         sys.exit(2)
 
 
--- a/pdb2pka/substruct/Algorithms.cpp
+++ b/pdb2pka/substruct/Algorithms.cpp
@@ -136,7 +136,7 @@
 
   for (int i=0; i<cliq.size; i++) {
     cerr << cliq.vertex[i] << " ";
-    PyList_Append(erg, PyInt_FromLong(cliq.vertex[i]));
+    PyList_Append(erg, PyLong_FromLong(cliq.vertex[i]));
   }
   
   for(int i = 0; i < n; i++) {
@@ -242,7 +242,7 @@
     PyObject* cli = PyList_New(0);
 
     for (int j=0; j < collector.cliques[i].size; j++) {
-      PyList_Append(cli, PyInt_FromLong(collector.cliques[i].vertex[j]));
+      PyList_Append(cli, PyLong_FromLong(collector.cliques[i].vertex[j]));
     }
 
     PyList_Append(result, cli);
@@ -299,10 +299,24 @@
 };
 
 /*----------------------------------------------------------------------------*/
+
+// Followed
+//   https://stackoverflow.com/questions/28305731/compiler-cant-find-py-initmodule-is-it-deprecated-and-if-so-what-should-i
+// to replace Py_InitModule
+
+static struct PyModuleDef cModMethods =
+{
+    PyModuleDef_HEAD_INIT,
+    MODULE_STR,  /* name of module */
+    "",          /* module documentation, may be NULL */
+    -1,          /* size of per-interpreter state of the module, or -1 if the module keeps state in global variables. */
+    ModuleMethods
+};
+
 PyMODINIT_FUNC MODULE_INIT()
 {
     
-    PyObject *ThisModule = Py_InitModule( MODULE_STR , ModuleMethods);
+    PyObject *ThisModule = PyModule_Create(&cModMethods);
     PyObject *ModuleDict = PyModule_GetDict(ThisModule);  
 
     // for correct handling of Numpy