diff --git a/propka/coupled_groups.py b/propka/coupled_groups.py
index aeed93e..bb47dfe 100644
--- a/propka/coupled_groups.py
+++ b/propka/coupled_groups.py
@@ -369,16 +369,17 @@ class NonCovalentlyCoupledGroups:
             "(difference: {diff_pka:>6.2f})\n"
             " Swapped pKa's:       {swap1:>6.2f}, {swap2:>6.2f}  "
             "(difference: {shift1:>6.2f}, {shift2:>6.2f})"
- ).format(
-     label1=group1.label, label2=group2.label,
-     coupl_fact=data['coupling_factor'], def_energy=data['default_energy'],
-     swap_energy=data['swapped_energy'],
-     diff_energy=data['default_energy']-data['swapped_energy'], ph=data['pH'],
-     int_energy=data['interaction_energy'], pka1=group1.intrinsic_pka,
-     pka2=group2.intrinsic_pka,
-     diff_pka=group1.intrinsic_pka-group2.intrinsic_pka, 
-     swap1=data['swapped_pka1'], swap2=data['swapped_pka2'],
-     shift1=data['pka_shift1'], shift2=data['pka_shift2'])
+        ).format(
+            label1=group1.label, label2=group2.label,
+            coupl_fact=data['coupling_factor'],
+            def_energy=data['default_energy'],
+            swap_energy=data['swapped_energy'],
+            diff_energy=data['default_energy']-data['swapped_energy'],
+            ph=data['pH'], int_energy=data['interaction_energy'],
+            pka1=group1.intrinsic_pka, pka2=group2.intrinsic_pka,
+            diff_pka=group1.intrinsic_pka-group2.intrinsic_pka,
+            swap1=data['swapped_pka1'], swap2=data['swapped_pka2'],
+            shift1=data['pka_shift1'], shift2=data['pka_shift2'])
         return str_
 
 
diff --git a/propka/group.py b/propka/group.py
index 125b0dc..9251275 100644
--- a/propka/group.py
+++ b/propka/group.py
@@ -473,11 +473,13 @@ class Group:
                     str_ += '*'
                 else:
                     str_ += ' '
-                str_ += " {0:4d}{1:>2s} ".format(int(100.0*self.buried), "%")
-                str_ += " %6.2lf %4d" % (self.energy_volume, self.num_volume)
-                str_ += " %6.2lf %4d" % (self.energy_local, self.num_local)
+                str_ += " {0:4d}{1:>2s} ".format(int(100.0*self.buried), "%%")
+                str_ += " {0:6.2f} {1:4d}".format(
+                    self.energy_volume, int(self.num_volume))
+                str_ += " {0:6.2f} {1:4d}".format(
+                    self.energy_local, int(self.num_local))
             else:
-                str_ += "%40s" % (" ")
+                str_ += "{0:>40s}".format(" ")
             # add the determinants
             for type_ in ['sidechain', 'backbone', 'coulomb']:
                 str_ += self.get_determinant_for_string(type_, line_number)
@@ -496,11 +498,11 @@ class Group:
             string
         """
         if number >= len(self.determinants[type_]):
-            empty_determinant = "%s%4d%2s" % ("XXX", 0, "X")
-            return "%8.2lf %s" % (0.0, empty_determinant)
+            return "    0.00 XXX   0 X"
         else:
             determinant = self.determinants[type_][number]
-            return "%8.2lf %s" % (determinant.value, determinant.label)
+            return "{0:8.2f} {1:s}".format(
+                determinant.value, determinant.label)
 
     def calculate_total_pka(self):
         """Calculate total pKa based on determinants associated with this
@@ -551,15 +553,16 @@ class Group:
         penalty = ''
         if self.coupled_titrating_group:
             penalty = (
-                ' NB: Discarded due to coupling with %s'
-                % self.coupled_titrating_group.label)
+                ' NB: Discarded due to coupling with {0:s}'.format(
+                    self.coupled_titrating_group.label))
         fmt = (
             "   {g.label:>9s} {g.pka_value:8.2f} {g.model_pka:10.2f} "
             "{type:>18s}   {penalty:s}\n")
         return fmt.format(g=self, type=ligand_type, penalty=penalty)
 
     def __str__(self):
-        return 'Group (%s) for %s' % (self.type, self.atom)
+        str_ = 'Group ({0:s}) for {1:s}'.format(self.type, str(self.atom))
+        return str_
 
     def calculate_folding_energy(self, parameters, ph=None, reference=None):
         """Return the electrostatic energy of this residue at specified pH.
@@ -1243,8 +1246,9 @@ def is_group(parameters, atom):
     elif parameters.ligand_typing == 'groups':
         ligand_group = is_ligand_group_by_groups(parameters, atom)
     else:
-        raise Exception('Unknown ligand typing method \'%s\''
-                        % parameters.ligand_typing)
+        raise Exception(
+            'Unknown ligand typing method \'{0.s}\''.format(
+                parameters.ligand_typing))
     if ligand_group:
         return ligand_group
     return None
@@ -1276,9 +1280,9 @@ def is_protein_group(parameters, atom):
         if atom.count_bonded_elements('O') == 1:
             return BBCGroup(atom)
     ### Filters for side chains based on PDB protein atom names
-    key = '%s-%s' % (atom.res_name, atom.name)
+    key = '{0:s}-{1:s}'.format(atom.res_name, atom.name)
     if key in parameters.protein_group_mapping.keys():
-        class_str = "%sGroup" % parameters.protein_group_mapping[key]
+        class_str = "{0:s}Group".format(parameters.protein_group_mapping[key])
         group_class = globals()[class_str]
         return group_class(atom)
     return None
diff --git a/propka/hybrid36.py b/propka/hybrid36.py
index 8acf0c1..b1e934a 100644
--- a/propka/hybrid36.py
+++ b/propka/hybrid36.py
@@ -20,7 +20,7 @@ def decode(input_string):
     Returns:
         integer
     """
-    value_error_message = "invalid literal for hybrid-36 conversion: '%s'"
+    value_error_message = "invalid literal for hybrid-36 conversion: '{0:s}'"
 
     original_input_string = input_string
     input_string = input_string.strip()
@@ -33,7 +33,7 @@ def decode(input_string):
         sign = 1
 
     if len(input_string) == 0:
-        raise ValueError(value_error_message % input_string)
+        raise ValueError(value_error_message.format(input_string))
 
     # See http://cci.lbl.gov/hybrid_36/ for documentation on the format.
 
@@ -49,13 +49,13 @@ def decode(input_string):
         reference = (16 * 36 ** (num_chars - 1) + 10 ** num_chars)
         _hybrid36_set = _HYBRID36_LOWER_SET
     else:
-        raise ValueError(value_error_message % original_input_string)
+        raise ValueError(value_error_message.format(original_input_string))
 
     # Check the validity of the input string: ASCII characters should be
     # either all uppercase or all lowercase.
     for char in input_string[1:]:
         if char not in _hybrid36_set:
-            raise ValueError(value_error_message % original_input_string)
+            raise ValueError(value_error_message.format(original_input_string))
 
     # Convert with the int function.
     return sign * (int(input_string, 36) + reference)
diff --git a/propka/iterative.py b/propka/iterative.py
index 5bafc8e..c3f999a 100644
--- a/propka/iterative.py
+++ b/propka/iterative.py
@@ -211,8 +211,9 @@ def add_determinants(iterative_interactions, version, _=None):
                 iteratives.append(new_iterative)
                 done_group.append(group)
     # Initialize iterative scheme
-    debug("\n   --- pKa iterations (%d groups, %d interactions) ---" %
-          (len(iteratives), len(iterative_interactions)))
+    debug(
+        "\n   --- pKa iterations ({0:d} groups, {1:d} interactions) ---".format(
+            len(iteratives), len(iterative_interactions)))
     converged = False
     iteration = 0
     # set non-iterative pka values as first step
@@ -262,18 +263,18 @@ def add_determinants(iterative_interactions, version, _=None):
             itres.pka_iter.append(itres.pka_new)
 
         if iteration == 10:
-            info("did not converge in %d iterations" % (iteration))
+            info("did not converge in {0:d} iterations".format(iteration))
             break
     # printing pKa iterations
     # formerly was conditioned on if options.verbosity >= 2 - now unnecessary
-    str_ = "%12s" % (" ")
+    str_ = '            '
     for index in range(iteration+1):
-        str_ += "%8d" % (index)
+        str_ += "{0:>8d}".format(index)
     debug(str_)
     for itres in iteratives:
-        str_ = "%s   " % (itres.label)
+        str_ = "{0:s}   ".format(itres.label)
         for pka in itres.pka_iter:
-            str_ += "%8.2lf" % (pka)
+            str_ += "{0:>8.2f}".format(pka)
         if not itres.converged:
             str_ += " *"
         debug(str_)
diff --git a/propka/lib.py b/propka/lib.py
index ab54850..6e45972 100644
--- a/propka/lib.py
+++ b/propka/lib.py
@@ -29,7 +29,7 @@ def open_file_for_reading(input_file):
     try:
         file_ = open(input_file, 'rt')
     except:
-        raise IOError('Cannot find file %s' % input_file)
+        raise IOError('Cannot find file {0:s}'.format(input_file))
     return file_
 
 
@@ -52,7 +52,7 @@ def open_file_for_writing(input_file):
     try:
         file_ = open(input_file, 'wt')
     except FileNotFoundError:
-        raise Exception('Could not open %s' % input_file)
+        raise Exception('Could not open {0:s}'.format(input_file))
     return file_
 
 
@@ -228,7 +228,7 @@ def build_parser(parser=None):
     group.add_argument(
         "-p", "--parameters", dest="parameters",
         default=pkg_resources.resource_filename(__name__, "propka.cfg"),
-        help="set the parameter file [%(default)s]")
+        help="set the parameter file [{default:s}]")
     try:
         group.add_argument(
             "--log-level",
@@ -306,7 +306,8 @@ def loadOptions(args):
             try:
                 chain, resnum, inscode = parse_res_string(res_str)
             except ValueError:
-                _LOGGER.critical('Invalid residue string: "%s"', res_str)
+                _LOGGER.critical(
+                    'Invalid residue string: "{0:s}"'.format(res_str))
                 sys.exit(1)
             res_list.append((chain, resnum, inscode))
         options.titrate_only = res_list
@@ -332,9 +333,9 @@ def make_tidy_atom_label(name, element):
         label = name
     else: # if less than 4 characters long, insert white space as needed
         if len(element) == 1:
-            label = ' %-3s' % name
+            label = ' {0:<3s}'.format(name)
         else: # The element should occupy the two first chars
-            label = '%-4s' % name
+            label = '{0:<4s}'.format(name)
     return label
 
 
@@ -365,7 +366,7 @@ def write_file(filename, lines):
     """
     file_ = open_file_for_writing(filename)
     for line in lines:
-        file_.write("%s\n" % (line))
+        file_.write("{0:s}\n".format(line))
     file_.close()
 
 
diff --git a/propka/ligand_pka_values.py b/propka/ligand_pka_values.py
index 43e9dac..d7a5423 100644
--- a/propka/ligand_pka_values.py
+++ b/propka/ligand_pka_values.py
@@ -42,7 +42,7 @@ class LigandPkaValues:
                               map(lambda dir: os.path.join(dir, program),
                                   path))]
         if len(locs) == 0:
-            str_ = "'Error: Could not find %s." % program
+            str_ = "'Error: Could not find {0:s}.".format(program)
             str_ += ' Please make sure that it is found in the path.'
             info(str_)
             sys.exit(-1)
@@ -73,7 +73,7 @@ class LigandPkaValues:
             max_ph:  maximum pH value
         """
         for name in molecule.conformation_names:
-            filename = '%s_%s' % (molecule.name, name)
+            filename = '{0:s}_{1:s}'.format(molecule.name, name)
             self.get_marvin_pkas_for_conformation_container(
                 molecule.conformations[name], name=filename,
                 reuse=molecule.options.reuse_ligand_mol2_file,
@@ -113,7 +113,7 @@ class LigandPkaValues:
         # do one molecule at the time so we don't confuse marvin
         molecules = propka.lib.split_atoms_into_molecules(atoms)
         for i, molecule in enumerate(molecules):
-            filename = '%s_%d.mol2' % (name, i+1)
+            filename = '{0:s}_{1:d}.mol2'.format(name, i+1)
             self.get_marvin_pkas_for_molecule(
                 molecule, filename=filename, reuse=reuse, num_pkas=num_pkas,
                 min_ph=min_ph, max_ph=max_ph)
@@ -136,8 +136,10 @@ class LigandPkaValues:
             propka.pdb.write_mol2_for_atoms(atoms, filename)
         # check that we actually have a file to work with
         if not os.path.isfile(filename):
-            errstr = ("Didn't find a user-modified file '%s' - generating one"
-                      % filename)
+            errstr = (
+                "Didn't find a user-modified file '{0:s}' "
+                "- generating one".format(
+                    filename))
             warning(errstr)
             propka.pdb.write_mol2_for_atoms(atoms, filename)
         # Marvin calculate pKa values
@@ -155,11 +157,11 @@ class LigandPkaValues:
                  '*********************************************')
             info('* Warning: Marvin execution failed:                        '
                  '                                            *')
-            info('* %-100s *' % errors)
+            info('* {0:<100s} *'.format(errors))
             info('*                                                          '
                  '                                            *')
             info('* Please edit the ligand mol2 file and re-run PropKa with '
-                 'the -l option: %29s *' % filename)
+                 'the -l option: {0:>29s} *'.format(filename))
             info('***********************************************************'
                  '*********************************************')
             sys.exit(-1)
@@ -169,7 +171,7 @@ class LigandPkaValues:
         for i, index in enumerate(indices):
             atoms[index].marvin_pka = pkas[i]
             atoms[index].charge = {'a': -1, 'b': 1}[types[i]]
-            info('%s model pKa: %.2f' % (atoms[index], pkas[i]))
+            info('{0:s} model pKa: {1:<.2f}'.format(atoms[index], pkas[i]))
 
     @staticmethod
     def extract_pkas(output):
@@ -188,7 +190,7 @@ class LigandPkaValues:
         values = values.split('\t')
         # format values
         types = [
-            tags[i][0] for i in range(1, len(tags)-1) 
+            tags[i][0] for i in range(1, len(tags)-1)
             if len(values) > i and values[i] != '']
         indices = [int(a)-1 for a in values[-1].split(',') if a != '']
         values = [float(v.replace(',', '.')) for v in values[1:-1] if v != '']
diff --git a/propka/molecular_container.py b/propka/molecular_container.py
index 0c97448..3d2409f 100644
--- a/propka/molecular_container.py
+++ b/propka/molecular_container.py
@@ -50,7 +50,8 @@ class Molecular_container:
             self.version = version_class(parameters)
         except AttributeError as err:
             print(err)
-            errstr = 'Error: Version %s does not exist' % parameters.version
+            errstr = 'Error: Version {0:s} does not exist'.format(
+                parameters.version)
             raise Exception(errstr)
         # read the input file
         if input_file_extension[0:4] == '.pdb':
@@ -89,7 +90,7 @@ class Molecular_container:
             # do some additional set up
             self.additional_setup_when_reading_input_file()
         else:
-            info('Unrecognized input file:%s' % input_file)
+            info('Unrecognized input file:{0:s}'.format(input_file))
             sys.exit(-1)
 
     def top_up_conformations(self):
@@ -152,8 +153,10 @@ class Molecular_container:
                 if group_to_add:
                     avr_group += group_to_add
                 else:
-                    str_ = ('Group %s could not be found in conformation %s.'
-                            % (group.atom.residue_label, name))
+                    str_ = (
+                        'Group {0:s} could not be found in '
+                        'conformation {0:s}.'.format(
+                            group.atom.residue_label, name))
                     warning(str_)
             # ... and store the average value
             avr_group = avr_group / len(self.conformation_names)
@@ -178,16 +181,16 @@ class Molecular_container:
             options:  options object
         """
         # write out the average conformation
-        filename = os.path.join('%s.pka' % (self.name))
+        filename = os.path.join('{0:s}.pka'.format(self.name))
         # if the display_coupled_residues option is true, write the results out
         # to an alternative pka file
         if self.options.display_coupled_residues:
-            filename = os.path.join('%s_alt_state.pka' % (self.name))
+            filename = os.path.join('{0:s}_alt_state.pka'.format(self.name))
         if (hasattr(self.version.parameters, 'output_file_tag')
                 and len(self.version.parameters.output_file_tag) > 0):
             filename = os.path.join(
-                '%s_%s.pka' % (self.name,
-                               self.version.parameters.output_file_tag))
+                '{0:s}_{1:s}.pka'.format(
+                    self.name, self.version.parameters.output_file_tag))
         propka.output.write_pka(
             self, self.version.parameters, filename=filename,
             conformation='AVR', reference=reference)
diff --git a/propka/output.py b/propka/output.py
index d2ed804..f533316 100644
--- a/propka/output.py
+++ b/propka/output.py
@@ -5,9 +5,9 @@ from propka.lib import info
 
 def print_header():
     """Print header section of output."""
-    str_ = "%s\n" % get_propka_header()
-    str_ += "%s\n" % get_references_header()
-    str_ += "%s\n" % get_warning_header()
+    str_ = "{0:s}\n".format(get_propka_header())
+    str_ += "{0:s}\n".format(get_references_header())
+    str_ += "{0:s}\n".format(get_warning_header())
     info(str_)
 
 
@@ -25,10 +25,10 @@ def write_pdb(protein, pdbfile=None, filename=None, include_hydrogens=False,
     if pdbfile is None:
         # opening file if not given
         if filename is None:
-            filename = "%s.pdb" % (protein.name)
+            filename = "{0:s}.pdb".format(protein.name)
         # TODO - this would be better as a context manager
         pdbfile = open(filename, 'w')
-        info("writing pdbfile %s" % (filename))
+        info("writing pdbfile {0:s}".format(filename))
         close_file = True
     else:
         # don't close the file, it was opened in a different place
@@ -67,20 +67,20 @@ def write_pka(protein, parameters, filename=None, conformation='1A',
     # TODO - the code immediately overrides the verbose argument; why?
     verbose = True
     if filename is None:
-        filename = "%s.pka" % (protein.name)
+        filename = "{0:s}.pka".format(protein.name)
     # TODO - this would be much better with a context manager
     file_ = open(filename, 'w')
     if verbose:
-        info("Writing %s" % (filename))
+        info("Writing {0:s}".format(filename))
     # writing propka header
-    str_ = "%s\n" % get_propka_header()
-    str_ += "%s\n" % get_references_header()
-    str_ += "%s\n" % get_warning_header()
+    str_ = "{0:s}\n".format(get_propka_header())
+    str_ += "{0:s}\n".format(get_references_header())
+    str_ += "{0:s}\n".format(get_warning_header())
     # writing pKa determinant section
     str_ += get_determinant_section(protein, conformation, parameters)
     # writing pKa summary section
     str_ += get_summary_section(protein, conformation, parameters)
-    str_ += "%s\n" % get_the_line()
+    str_ += "{0:s}\n".format(get_the_line())
     # printing Folding Profile
     str_ += get_folding_profile_section(
         protein, conformation=conformation, reference=reference,
@@ -109,17 +109,18 @@ def print_tm_profile(protein, reference="neutral", window=[0., 14., 1.],
         _:  Boolean for verbosity
         options:  options object
     """
-    profile = protein.getTmProfile(reference=reference, grid=[0., 14., 0.1],
-                                   tms=tms, ref=ref, options=options)
+    profile = protein.getTmProfile(
+        reference=reference, grid=[0., 14., 0.1], tms=tms, ref=ref,
+        options=options)
     if profile is None:
         str_ = "Could not determine Tm-profile\n"
     else:
-        str_ = " suggested Tm-profile for %s\n" % (protein.name)
+        str_ = " suggested Tm-profile for {0:s}\n".format(protein.name)
         for (ph, tm_) in profile:
             if (ph >= window[0] and ph <= window[1]
-                    and (ph%window[2] < 0.01
-                         or ph%window[2] > 0.99*window[2])):
-                str_ += "%6.2lf%10.2lf\n" % (ph, tm_)
+                    and (ph % window[2] < 0.01
+                         or ph % window[2] > 0.99*window[2])):
+                str_ += "{0:>6.2f}{1:>10.2f}\n".format(ph, tm_)
         info(str_)
 
 
@@ -160,7 +161,7 @@ def get_determinant_section(protein, conformation, parameters):
         string
     """
     # getting the same order as in propka2.0
-    str_ = "%s\n" % get_determinants_header()
+    str_ = "{0:s}\n".format(get_determinants_header())
     # printing determinants
     for chain in protein.conformations[conformation].chains:
         for residue_type in parameters.write_out_order:
@@ -169,8 +170,9 @@ def get_determinant_section(protein, conformation, parameters):
                 if g.atom.chain_id == chain]
             for group in groups:
                 if group.residue_type == residue_type:
-                    str_ += "%s" % group.get_determinant_string(
-                        parameters.remove_penalised_group)
+                    str_ += "{0:s}".format(
+                        group.get_determinant_string(
+                            parameters.remove_penalised_group))
     # Add a warning in case of coupled residues
     if (protein.conformations[conformation].non_covalently_coupled_groups
             and not protein.options.display_coupled_residues):
@@ -190,13 +192,14 @@ def get_summary_section(protein, conformation, parameters):
     Returns:
         string
     """
-    str_ = "%s\n" % get_summary_header()
+    str_ = "{0:s}\n".format(get_summary_header())
     # printing pKa summary
     for residue_type in parameters.write_out_order:
         for group in protein.conformations[conformation].groups:
             if group.residue_type == residue_type:
-                str_ += "%s" % group.get_summary_string(
-                    parameters.remove_penalised_group)
+                str_ += "{0:s}".format(
+                    group.get_summary_string(
+                        parameters.remove_penalised_group))
     return str_
 
 
@@ -219,8 +222,8 @@ def get_folding_profile_section(protein, conformation='AVR',
     """
     str_ = get_the_line()
     str_ += "\n"
-    str_ += "Free energy of %9s (kcal/mol) as a function" % direction
-    str_ += " of pH (using %s reference)\n" %  reference
+    str_ += "Free energy of {0:>9s} (kcal/mol) as a function".format(direction)
+    str_ += " of pH (using {0:s} reference)\n".format(reference)
     profile, [ph_opt, dg_opt], [dg_min, dg_max], [ph_min, ph_max] = (
         protein.get_folding_profile(
             conformation=conformation, reference=reference,
@@ -230,27 +233,28 @@ def get_folding_profile_section(protein, conformation='AVR',
     else:
         for (ph, dg) in profile:
             if ph >= window[0] and ph <= window[1]:
-                if ph%window[2] < 0.05 or ph%window[2] > 0.95:
-                    str_ += "%6.2lf%10.2lf\n" % (ph, dg)
+                if ph % window[2] < 0.05 or ph % window[2] > 0.95:
+                    str_ += "{0:>6.2f}{1:>10.2f}\n".format(ph, dg)
         str_ += "\n"
     if ph_opt is None or dg_opt is None:
         str_ += "Could not determine pH optimum\n"
     else:
-        str_ += "The pH of optimum stability is %4.1lf" % ph_opt
-        str_ += (" for which the free energy is %6.1lf kcal/mol at 298K\n"
-                 % dg_opt)
+        str_ += "The pH of optimum stability is {0:>4.1f}".format(ph_opt)
+        str_ += (
+            " for which the free energy is {0:>6.1f} kcal/mol at 298K\n".format(
+                dg_opt))
     if dg_min is None or dg_max is None:
         str_ += "Could not determine pH values where the free energy"
-        str_ += " is within 80 %s of minimum\n" % ("%")
+        str_ += " is within 80 %% of minimum\n"
     else:
         str_ += "The free energy is within 80 %% of maximum"
-        str_ += " at pH %4.1lf to %4.1lf\n" % (dg_min, dg_max)
+        str_ += " at pH {0:>4.1f} to {1:>4.1f}\n".format(dg_min, dg_max)
     if ph_min is None or ph_max is None:
         str_ += "Could not determine the pH-range where the free"
         str_ += " energy is negative\n\n"
     else:
         str_ += "The free energy is negative in the range"
-        str_ += " %4.1lf - %4.1lf\n\n" % (ph_min, ph_max)
+        str_ += " {0:>4.1f} - {1:>4.1f}\n\n".format(ph_min, ph_max)
     return str_
 
 
@@ -278,7 +282,7 @@ def get_charge_profile_section(protein, conformation='AVR', _=None):
     if pi_pro is None or pi_mod is None:
         str_ += "Could not determine the pI\n\n"
     else:
-        str_ += ("The pI is %5.2lf (folded) and %5.2lf (unfolded)\n")
+        str_ += ("The pI is {0:>5.2f} (folded) and {1:>5.2f} (unfolded)\n")
     return str_
 
 
@@ -290,7 +294,8 @@ def write_jackal_scap_file(mutation_data=None, filename="1xxx_scap.list",
     """
     with open(filename, 'w') as file_:
         for chain_id, _, res_num, code2 in mutation_data:
-            str_ = "%s, %d, %s\n" % (chain_id, res_num, code2)
+            str_ = "{chain:s}, {num:d}, {code:s}\n".format(
+                chain=chain_id, num=res_num, code=code2)
         file_.write(str_)
 
 
@@ -302,9 +307,9 @@ def write_scwrl_sequence_file(sequence, filename="x-ray.seq", _=None):
     with open(filename, 'w') as file_:
         start = 0
         while len(sequence[start:]) > 60:
-            file_.write("%s\n" % (sequence[start:start+60]))
+            file_.write("{0:s}s\n".format(sequence[start:start+60]))
             start += 60
-        file_.write("%s\n" % (sequence[start:]))
+        file_.write("{0:s}\n".format(sequence[start:]))
 
 
 def get_propka_header():
@@ -314,7 +319,7 @@ def get_propka_header():
         string
     """
     today = date.today()
-    str_ = "propka3.1 %93s\n" % (today)
+    str_ = "propka3.1 {0:>93s}\n".format(today)
     str_ += ("---------------------------------------------------------------"
              "----------------------------------------\n")
     str_ += ("--                                                             "
@@ -468,19 +473,19 @@ def make_interaction_map(name, list_, interaction):
         for i, group1 in enumerate(list_):
             for group2 in list_[i:]:
                 if interaction(group1, group2):
-                    res += 'Coupling: %9s - %9s\n' % (group1.label,
-                                                      group2.label)
+                    res += 'Coupling: {0:>9s} - {1:>9s}\n'.format(
+                        group1.label, group2.label)
         return res
     # Name and map header
-    res = '%s\n%12s' % (name, '')
+    res = '{0:s}\n{1:>12s}'.format(name, '')
     for group in list_:
-        res += '%9s | ' % group.label
+        res += '{0:>9s} | '.format(group.label)
     # do the map
     for group1 in list_:
-        res += '\n%-12s' % (group1.label)
+        res += '\n{0:<12s}'.format(group1.label)
         for group2 in list_:
             tag = ''
             if interaction(group1, group2):
                 tag = '    X     '
-            res += '%10s| '%tag
+            res += '{0:>10s}| '.format(tag)
     return res
diff --git a/propka/parameters.py b/propka/parameters.py
index 4bdff07..960cb81 100644
--- a/propka/parameters.py
+++ b/propka/parameters.py
@@ -163,7 +163,8 @@ class Parameters:
         value = float(words[1])
         setattr(self, words[0], value)
         value_sq = value*value
-        setattr(self, "%s_squared" % words[0], value_sq)
+        attr = "{0:s}_squared".format(words[0])
+        setattr(self, attr, value_sq)
 
     def parse_parameter(self, words):
         """Parse field to parameters.
@@ -206,7 +207,7 @@ class Parameters:
         """Print interaction parameters."""
         info('--------------- Model pKa values ----------------------')
         for k in self.model_pkas:
-            info('%3s %8.2f' % (k, self.model_pkas[k]))
+            info('{0:>3s} {1:8.2f}'.format(k, self.model_pkas[k]))
 
         info('')
         info('--------------- Interactions --------------------------')
@@ -227,10 +228,10 @@ class Parameters:
             for group2 in lgroups:
                 fmt = "{grp1:>3s} {grp2:>3s} {mat:1s} {val1:4} {val2:4}"
                 interaction = fmt.format(
-                        grp1=group1, grp2=group2,
-                        mat=self.interaction_matrix[group1][group2],
-                        val1=self.sidechain_cutoffs.get_value(group1, group2)[0],
-                        val2=self.sidechain_cutoffs.get_value(group1, group2)[1])
+                    grp1=group1, grp2=group2,
+                    mat=self.interaction_matrix[group1][group2],
+                    val1=self.sidechain_cutoffs.get_value(group1, group2)[0],
+                    val2=self.sidechain_cutoffs.get_value(group1, group2)[1])
                 map_interaction = ''
                 if group2 in map_:
                     for val in map_[group2]:
@@ -243,21 +244,21 @@ class Parameters:
                                 != self.interaction_matrix[group1][group2]):
                             map_interaction += '* '
                         if (self.sidechain_cutoffs.get_value(group1, val)[0]
-                                != self.sidechain_cutoffs.get_value(group1,
-                                                                    group2)[0]
-                                or self.sidechain_cutoffs.get_value(group1,
-                                                                    val)[1]
-                                != self.sidechain_cutoffs.get_value(group1,
-                                                                    group2)[1]):
+                                != self.sidechain_cutoffs.get_value(
+                                    group1, group2)[0]
+                                or self.sidechain_cutoffs.get_value(
+                                    group1, val)[1]
+                                != self.sidechain_cutoffs.get_value(
+                                    group1, group2)[1]):
                             map_interaction += '! '
                         else:
                             map_interaction += '  '
                     if (len(map_[group2]) == 0
-                            and (self.sidechain_cutoffs.get_value(group1,
-                                                                  group2)[0]
+                            and (self.sidechain_cutoffs.get_value(
+                                group1, group2)[0]
                                  != 3
-                                 or self.sidechain_cutoffs.get_value(group1,
-                                                                     group2)[1]
+                                 or self.sidechain_cutoffs.get_value(
+                                     group1, group2)[1]
                                  != 4)):
                         map_interaction += '?  '
                 info(interaction, map_interaction)
@@ -306,30 +307,32 @@ O2
         lgroups = ['CG', 'C2N', 'N30', 'N31', 'N32', 'N33', 'NAR', 'OCO',
                    'NP1', 'OH', 'O3', 'CL', 'F', 'NAM', 'N1', 'O2', 'OP',
                    'SH']
-        str_ = """
-\\begin{longtable}{lllll}
-\\caption{Ligand interaction parameters. For interactions not listed, the default value of %s is applied.}
-\\label{tab:ligand_interaction_parameters}\\\\
-
-\\toprule
-Group1 & Group2 & Interaction & c1 &c2 \\\\
-\\midrule
-\\endfirsthead
-
-\\multicolumn{5}{l}{\\emph{continued from the previous page}}\\\\
-\\toprule
-Group1 & Group2 & Interaction & c1 &c2 \\\\
-\\midrule
-\\endhead
-
-\\midrule
-\\multicolumn{5}{r}{\\emph{continued on the next page}}\\\\
-\\endfoot
-
-\\bottomrule
-\\endlastfoot
-
-""" % (self.sidechain_cutoffs.default)
+        lines = [
+            "",
+            "\\begin{{longtable}}{{lllll}}",
+            ("\\caption{{Ligand interaction parameters. For interactions not "
+             "listed, the default value of {0:s} is applied.}}").format(
+                 self.sidechain_cutoffs.default),
+            "\\label{{tab:ligand_interaction_parameters}}\\\\",
+            "\\toprule",
+            "Group1 & Group2 & Interaction & c1 &c2 \\\\",
+            "\\midrule",
+            "\\endfirsthead",
+            "",
+            "\\multicolumn{{5}}{{l}}{\\emph{{continued from the previous page}}}\\\\",
+            "\\toprule",
+            "Group1 & Group2 & Interaction & c1 &c2 \\\\",
+            "\\midrule",
+            "\\endhead",
+            "",
+            "\\midrule",
+            "\\multicolumn{{5}}{{r}}{\\emph{{continued on the next page}}}\\\\",
+            "\\endfoot",
+            "",
+            "\\bottomrule",
+            "\\endlastfoot",
+            ""]
+        str_ = "\n".join(lines)
         for group1 in agroups:
             for group2 in lgroups:
                 if self.interaction_matrix[group1][group2] == '-':
@@ -347,7 +350,7 @@ Group1 & Group2 & Interaction & c1 &c2 \\\\
                     self.sidechain_cutoffs.get_value(group1, group2)[1])
                 if group1 == group2:
                     break
-        str_ += '  \\end{longtable}\n'
+        str_ += '  \\end{{longtable}}\n'
         info(str_)
 
     def print_interactions_latex(self):
@@ -357,42 +360,43 @@ Group1 & Group2 & Interaction & c1 &c2 \\\\
                    'ARG', 'TRP', 'ROH', 'CG', 'C2N', 'N30', 'N31', 'N32',
                    'N33', 'NAR', 'OCO', 'NP1', 'OH', 'O3', 'CL', 'F', 'NAM',
                    'N1', 'O2', 'OP', 'SH']
-        str_ = """
-\\begin{longtable}{%s}
-\\caption{Ligand interaction parameters. For interactions not listed, the default value of %s is applied.}
-\\label{tab:ligand_interaction_parameters}\\\\
-
-\\toprule
-Group1 & Group2 & Interaction & c1 &c2 \\\\
-\\midrule
-\\endfirsthead
-
-\\multicolumn{5}{l}{\\emph{continued from the previous page}}\\\\
-\\toprule
-Group1 & Group2 & Interaction & c1 &c2 \\\\
-\\midrule
-\\endhead
-
-\\midrule
-\\multicolumn{5}{r}{\\emph{continued on the next page}}\\\\
-\\endfoot
-
-\\bottomrule
-\\endlastfoot
-
-""" % ('l'*len(agroups), self.sidechain_cutoffs.default)
+        lines = [
+            "",
+            "\\begin{longtable}{{{0:s}}}".format('l'*len(agroups)),
+            ("\\caption{{Ligand interaction parameters. For interactions not "
+             "listed, the default value of {0:s} is applied.}}").format(
+                 str(self.sidechain_cutoffs.default)),
+            "\\label{{tab:ligand_interaction_parameters}}\\\\",
+            "\\toprule",
+            "Group1 & Group2 & Interaction & c1 &c2 \\\\",
+            "\\midrule",
+            "\\endfirsthead",
+            "",
+            "\\multicolumn{{5}}{{l}}{\\emph{{continued from the previous page}}}\\\\",
+            "\\toprule",
+            "Group1 & Group2 & Interaction & c1 &c2 \\\\",
+            "\\midrule",
+            "\\endhead",
+            "",
+            "\\midrule",
+            "\\multicolumn{{5}}{{r}}{\\emph{{continued on the next page}}}\\\\",
+            "\\endfoot",
+            "",
+            "\\bottomrule",
+            "\\endlastfoot",
+            ""
+        ]
+        str_ = "\n".join(lines)
         for group1 in agroups:
             for group2 in agroups:
-                str_ += ('%3s & %3s & %1s & %4s & %4s\\\\ \n'
-                         % (group1, group2,
-                            self.interaction_matrix[group1][group2],
-                            self.sidechain_cutoffs.get_value(
-                                group1, group2)[0],
-                            self.sidechain_cutoffs.get_value(
-                                group1, group2)[1]))
+                fmt = '{g1:>3s} & {g2:>3s} & {mat:1s} & {val1:>4s} & {val2:>4s}\\\\ \n'
+                str_ += fmt.format(
+                    group1, group2, self.interaction_matrix[group1][group2],
+                    str(self.sidechain_cutoffs.get_value(group1, group2)[0]),
+                    str(self.sidechain_cutoffs.get_value(group1, group2)[1]))
                 if group1 == group2:
                     break
-        str_ += '  \\end{longtable}\n'
+        str_ += '  \\end{{longtable}}\n'
         info(str_)
 
 
@@ -450,7 +454,8 @@ class InteractionMatrix:
             group:  group to get
         """
         if group not in self.dictionary.keys():
-            str_ = '%s not found in interaction matrix %s' % (group, self.name)
+            str_ = '{0:s} not found in interaction matrix {1:s}'.format(
+                group, self.name)
             raise KeyError(str_)
         return self.dictionary[group]
 
@@ -465,12 +470,12 @@ class InteractionMatrix:
     def __str__(self):
         str_ = '      '
         for key in self.ordered_keys:
-            str_ += '%3s ' % key
+            str_ += '{0:>3s} '.format(key)
         str_ += '\n'
         for key1 in self.ordered_keys:
-            str_ += '%3s ' % key1
+            str_ += '{0:>3s} '.format(key1)
             for key2 in self.ordered_keys:
-                str_ += '%3s ' % self[key1][key2]
+                str_ += '{0:>3s} '.format(self[key1][key2])
             str_ += '\n'
         return str_
 
@@ -517,8 +522,9 @@ class PairwiseMatrix:
         """
         if key1 in self.dictionary and key2 in self.dictionary[key1]:
             if key1 != key2:
-                str_ = ('Parameter value for %s, %s defined more than once'
-                        % (key1, key2))
+                str_ = (
+                    'Parameter value for {0:s}, {1:s} defined more '
+                    'than once'.format(key1, key2))
                 warning(str_)
         if not key1 in self.dictionary:
             self.dictionary[key1] = {}
@@ -547,7 +553,8 @@ class PairwiseMatrix:
             matrix information
         """
         if group not in self.dictionary.keys():
-            str_ = '%s not found in interaction matrix %s' % (group, self.name)
+            str_ = '{0:s} not found in interaction matrix {1:s}'.format(
+                group, self.name)
             raise KeyError(str_)
         return self.dictionary[group]
 
@@ -563,5 +570,6 @@ class PairwiseMatrix:
         str_ = ''
         for key1 in self.keys():
             for key2 in self[key1].keys():
-                str_ += '%s %s %s\n' % (key1, key2, self[key1][key2])
+                str_ += '{0:s} {1:s} {2:s}\n'.format(
+                    key1, key2, self[key1][key2])
         return str_
diff --git a/propka/pdb.py b/propka/pdb.py
index 7d33299..0074938 100644
--- a/propka/pdb.py
+++ b/propka/pdb.py
@@ -58,23 +58,25 @@ def protein_precheck(conformations, names):
                     atoms_by_residue[res_id] = [atom]
         for res_id, res_atoms in atoms_by_residue.items():
             res_name = res_atoms[0].res_name
-            residue_label = '%3s%5s'%(res_name, res_id)
+            residue_label = '{0:>3s}{1:>5s}'.format(res_name, res_id)
             # ignore ligand residues
             if res_name not in EXPECTED_ATOM_NUMBERS:
                 continue
             # check for c-terminal
             if 'C-' in [a.terminal for a in res_atoms]:
                 if len(res_atoms) != EXPECTED_ATOM_NUMBERS[res_name]+1:
-                    str_ = ("Unexpected number (%d) of atoms in residue %s "
-                            "in conformation %s" % (len(res_atoms),
-                                                    residue_label, name))
+                    str_ = ("Unexpected number ({num:d}) of atoms in residue "
+                            "{res:s} in conformation {conf:s}".format(
+                                num=len(res_atoms), res=residue_label,
+                                conf=name))
                     warning(str_)
                 continue
             # check number of atoms in residue
             if len(res_atoms) != EXPECTED_ATOM_NUMBERS[res_name]:
-                str_ = ('Unexpected number (%d) of atoms in residue %s '
-                        'in conformation %s' % (len(res_atoms),
-                                                residue_label, name))
+                str_ = ("Unexpected number ({num:d}) of atoms in residue "
+                        "{res:s} in conformation {conf:s}".format(
+                            num=len(res_atoms), res=residue_label,
+                            conf=name))
                 warning(str_)
 
 
@@ -86,7 +88,8 @@ def resid_from_atom(atom):
     Returns
         string
     """
-    return '%4d %s %s' % (atom.res_num, atom.chain_id, atom.icode)
+    return '{0:>4d} {1:s} {2:s}'.format(
+        atom.res_num, atom.chain_id, atom.icode)
 
 
 def get_atom_lines_from_pdb(pdb_file, ignore_residues=[], keep_protons=False,
@@ -136,7 +139,7 @@ def get_atom_lines_from_pdb(pdb_file, ignore_residues=[], keep_protons=False,
                 alt_conf_tag = chr(ord(alt_conf_tag)+16)
             if alt_conf_tag == ' ':
                 alt_conf_tag = 'A'
-            conformation = '%d%s'%(model, alt_conf_tag)
+            conformation = '{0:d}{1:s}'.format(model, alt_conf_tag)
             # set the terminal
             if  tag == 'ATOM  ':
                 if (residue_name.strip() == 'N'
@@ -190,7 +193,7 @@ def write_mol2_for_atoms(atoms, filename):
         filename:  name of file
     """
     # TODO - header needs to be converted to format string
-    header = '@<TRIPOS>MOLECULE\n\n%d %d\nSMALL\nUSER_CHARGES\n'
+    header = '@<TRIPOS>MOLECULE\n\n{natom:d} {id:d}\nSMALL\nUSER_CHARGES\n'
     atoms_section = '@<TRIPOS>ATOM\n'
     for i, atom in enumerate(atoms):
         atoms_section += atom.make_mol2_line(i+1)
@@ -200,15 +203,16 @@ def write_mol2_for_atoms(atoms, filename):
         for j, atom2 in enumerate(atoms, i+1):
             if atom1 in atom2.bonded_atoms:
                 type_ = get_bond_order(atom1, atom2)
-                bonds_section += '%7d %7d %7d %7s\n' % (id_, i+1, j+1, type_)
+                bonds_section += '{0:>7d} {1:>7d} {2:>7d} {3:>7s}\n'.format(
+                    id_, i+1, j+1, type_)
                 id_ += 1
     substructure_section = '@<TRIPOS>SUBSTRUCTURE\n\n'
     if len(atoms) > 0:
-        substructure_section = ('@<TRIPOS>SUBSTRUCTURE\n%-7d %10s %7d\n'
-                                % (atoms[0].res_num, atoms[0].res_name,
-                                   atoms[0].numb))
+        substructure_section = (
+            '@<TRIPOS>SUBSTRUCTURE\n{0:<7d} {1:>10s} {2:>7d}\n'.format(
+                atoms[0].res_num, atoms[0].res_name, atoms[0].numb))
     out = propka.lib.open_file_for_writing(filename)
-    out.write(header % (len(atoms), id_-1))
+    out.write(header.format(natom=len(atoms), id=id_-1))
     out.write(atoms_section)
     out.write(bonds_section)
     out.write(substructure_section)
@@ -232,7 +236,7 @@ def get_bond_order(atom1, atom2):
     if '.ar' in atom2.sybyl_type:
         pi_electrons2 -= 1
     if pi_electrons1 > 0 and pi_electrons2 > 0:
-        type_ = '%d' % (min(pi_electrons1, pi_electrons2)+1)
+        type_ = '{0:d}'.format(min(pi_electrons1, pi_electrons2)+1)
     if '.ar' in atom1.sybyl_type and '.ar' in atom2.sybyl_type:
         type_ = 'ar'
     return type_
@@ -247,7 +251,7 @@ def write_input(molecular_container, filename):
     """
     out = propka.lib.open_file_for_writing(filename)
     for conformation_name in molecular_container.conformation_names:
-        out.write('MODEL %s\n' % conformation_name)
+        out.write('MODEL {0:s}\n'.format(conformation_name))
         # write atoms
         for atom in molecular_container.conformations[conformation_name].atoms:
             out.write(atom.make_input_line())
diff --git a/propka/protonate.py b/propka/protonate.py
index 8218176..cd0f6c8 100644
--- a/propka/protonate.py
+++ b/propka/protonate.py
@@ -71,7 +71,7 @@ class Protonate:
         """
         # atom is a protein atom
         if atom.type == 'atom':
-            key = '%3s-%s' % (atom.res_name, atom.name)
+            key = '{0:3s}-{1:s}'.format(atom.res_name, atom.name)
             if atom.terminal:
                 debug(atom.terminal)
                 key = atom.terminal