diff --git a/MANIFEST.in b/MANIFEST.in
index c58a255b9..8aeafdcd4 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -9,6 +9,7 @@ include prody/*/*.cpp
 include prody/*/*/*.cpp
 include prody/*/*.h
 include prody/*/*/*.h
+include prody/*/*/*/*.so
 include prody/tests/*/*.py
 include prody/tests/datafiles/*.coo
 include prody/tests/datafiles/*.dcd
@@ -20,3 +21,4 @@ include prody/utilities/datafiles/*.dat
 include prody/utilities/datafiles/*.txt
 include scripts/prody*
 include scripts/evol*
+include prody/protein/tabulated_energies.txt
\ No newline at end of file
diff --git a/docs/about/people.rst b/docs/about/people.rst
index a46c0516a..a092a2098 100644
--- a/docs/about/people.rst
+++ b/docs/about/people.rst
@@ -21,10 +21,12 @@ Development Team
 `James Krieger`_ was helping develop *ProDy* from 2017 and became the main 
 overseer and developer in mid 2020.
 
-`Hongchun Li`_ is currently maintaining and developing ANM and GNM servers, 
+`Hongchun Li`_ has helped maintain and develop ANM and GNM servers, 
 and made significant contributions to :mod:`.database` and :mod:`.dynamics` 
 including *SignDy* and Adaptive ANM.
 
+`JiYoung Lee`_ is the main developer of :mod:`.Pharmmaker`, for constructing pharmacophore models using the outputs from :mod:`.DruGUI`.
+
 `Yan Zhang`_ contributed significantly to the development of 
 the *cryo-EM* module, :mod:`.protein.emdmap`.
 
@@ -32,9 +34,9 @@ the *cryo-EM* module, :mod:`.protein.emdmap`.
 :mod:`.domain_decomposition`, :mod:`.dynamics.essa`, and
 :mod:`.dynamics.clustenm`. 
 
-`Karolina Mikulska-Ruminska`_ contributed significantly to the development of 
+`Karolina Mikulska-Ruminska`_ is one of the main developers since 2022, with the addition of the new modules 
 :mod:`.protein.interactions` (*InSty*), :mod:`.protein.waterbridges`
-(*WatFinder*), and :mod:`.dynamics.mechstiff` (*MechStiff*).
+(*WatFinder*), in addition to being the developer of :mod:`.dynamics.mechstiff` (*MechStiff*).
 
 `Anthony Bogetti`_ is overseeing the overall development of *ProDy* since 
 2024.
@@ -50,6 +52,12 @@ Blocks and Membrane ENM.
 `Lidio Meireles`_ provided insightful comments on the design of *ProDy*,
 and contributed to the development of :ref:`prody-apps`.
 
+`Mustafa Tekpinar`_ contributed to the development of MechStiff.
+
+`Luca Ponzoni`_ contributed to the development of SignDy.
+
+`David Koes`_ contributed to the development of drug discovery tools.
+
 Contributors
 ------------
 
@@ -67,7 +75,6 @@ contributions and feedback from the following individuals:
 insights.
 
 
-
 .. _Ahmet Bakan: https://scholar.google.com/citations?user=-QAYVgMAAAAJ&hl=en
 .. _Cihan Kaya: https://www.linkedin.com/in/cihan-kaya/
 .. _Bahar Lab: http://www.bahargroup.org/Faculty/bahar/
@@ -75,9 +82,9 @@ insights.
 .. _Anindita Dutta: http://www.linkedin.com/pub/anindita-dutta/5a/568/a90
 .. _Wenzhi Mao: http://www.linkedin.com/pub/wenzhi-mao/2a/29a/29
 .. _Lidio Meireles: http://www.linkedin.com/in/lidio
-.. _Ying Liu: http://www.linkedin.com/pub/ying-liu/15/48b/5a9
+.. _Ying Liu: https://www.linkedin.com/in/yingliu03/
 .. _Kian Ho: https://github.com/kianho
-.. _Gökçen Eraslan: http://blog.yeredusuncedernegi.com/
+.. _Gökçen Eraslan: https://github.com/gokceneraslan
 .. _Tim Lezon: https://scholar.google.pl/citations?user=1MwNI3EAAAAJ&hl=pl&oi=ao
 .. _Chakra Chennubhotla: http://www.csb.pitt.edu/Faculty/Chakra/
 .. _She (John) Zhang: https://www.linkedin.com/in/she-zhang-49164399/
@@ -87,4 +94,8 @@ insights.
 .. _Burak Kaynak: https://scholar.google.pl/citations?user=gP8RokwAAAAJ&hl=pl&oi=ao
 .. _Karolina Mikulska-Ruminska: https://scholar.google.pl/citations?user=IpyPHRwAAAAJ&hl=pl
 .. _Anthony Bogetti: https://scholar.google.pl/citations?hl=pl&user=9qQClIcAAAAJ
-.. _Frane Doljanin: https://github.com/fdoljanin
\ No newline at end of file
+.. _Frane Doljanin: https://github.com/fdoljanin
+.. _JiYoung Lee: https://scholar.google.com/citations?user=odKQmZcAAAAJ&hl=en
+.. _David Koes: https://bits.csb.pitt.edu/
+.. _Luca Ponzoni: https://scholar.google.it/citations?user=8vfPOYUAAAAJ&hl=en
+.. _Mustafa Tekpinar: https://scholar.google.com/citations?user=qeVv6o8AAAAJ&hl=en   
diff --git a/docs/devel/website.rst b/docs/devel/website.rst
index e1ef59965..3de4d9bb4 100644
--- a/docs/devel/website.rst
+++ b/docs/devel/website.rst
@@ -14,7 +14,7 @@ This is a short guide for building the ProDy website.
 Environment Setup
 --------------
 
-First log in to the ProDy webserver (prody.csb.pitt.edu) then run the following::
+First log in to your ProDy webserver and then run the following::
 
   $ conda deactivate
 
@@ -34,7 +34,7 @@ ProDy-website-workdir. You can then copy files back over afterwards.
 
 It's recommended to have the symbolic link called test_prody pointing to 
 your build directory instead and then you can monitor changes by going to 
-http://prody.csb.pitt.edu/test_prody/_build/html/ in your web browser.
+http://yourdomainname/test_prody/_build/html/ in your web browser.
 
 
 Updating from GitHub
diff --git a/docs/docs b/docs/docs
index 945c9b46d..7099f6ab8 120000
--- a/docs/docs
+++ b/docs/docs
@@ -1 +1 @@
-.
\ No newline at end of file
+ProDy/docs
\ No newline at end of file
diff --git a/docs/reference/proteins/waterbridges.rst b/docs/reference/proteins/waterbridges.rst
index a41a8e09b..148fab777 100644
--- a/docs/reference/proteins/waterbridges.rst
+++ b/docs/reference/proteins/waterbridges.rst
@@ -3,4 +3,4 @@ Water bridge finder (WatFinder)
 
 .. automodule:: prody.proteins.waterbridges
    :members:
-   :undoc-members:
\ No newline at end of file
+   :undoc-members:
diff --git a/docs/release/index.rst b/docs/release/index.rst
index cfe0dfaff..f838d806b 100644
--- a/docs/release/index.rst
+++ b/docs/release/index.rst
@@ -9,7 +9,8 @@ Release Notes
 .. toctree::
    :maxdepth: 2
    :glob:
-   
+
+   v2.5_series
    v2.4_series
    v2.3_series
    v2.2_series
diff --git a/prody/chromatin/functions.py b/prody/chromatin/functions.py
index a9ed3a472..314fea2f3 100644
--- a/prody/chromatin/functions.py
+++ b/prody/chromatin/functions.py
@@ -212,7 +212,7 @@ def showEmbedding(modes, labels=None, trace=True, headtail=True, cmap='prism'):
         X, Y, Z = V[:,:3].T
         
         f = figure()
-        ax = f.add_subplot(projection="3d")
+        ax = f.add_subplot(1,1,1,projection="3d")
         if trace:
             ax.plot(X, Y, Z, ':', color=[0.3, 0.3, 0.3])
         if labels is None:
@@ -240,4 +240,4 @@ def getDomainList(labels):
     ends = sites[1:]
     domains = np.array([starts, ends]).T
 
-    return domains
\ No newline at end of file
+    return domains
diff --git a/prody/dynamics/plotting.py b/prody/dynamics/plotting.py
index 3a87ee4dc..867933572 100644
--- a/prody/dynamics/plotting.py
+++ b/prody/dynamics/plotting.py
@@ -110,7 +110,7 @@ def showEllipsoid(modes, onto=None, n_std=2, scale=1., *args, **kwargs):
             show = child
             break
     if show is None:
-        show = cf.add_subplot(projection="3d")
+        show = cf.add_subplot(111,projection="3d")
     show.plot_wireframe(x, y, z, rstride=6, cstride=6, *args, **kwargs)
     if onto is not None:
         onto = list(onto)
@@ -421,7 +421,7 @@ def showProjection(ensemble=None, modes=None, projection=None, *args, **kwargs):
                 show = child
                 break
         if show is None:
-            show = cf.add_subplot(projection="3d")
+            show = cf.add_subplot(111,projection="3d")
         plot = show.plot
         text = show.text
 
diff --git a/prody/dynamics/sampling.py b/prody/dynamics/sampling.py
index dd2c66453..272b230bd 100644
--- a/prody/dynamics/sampling.py
+++ b/prody/dynamics/sampling.py
@@ -51,9 +51,10 @@ def sampleModes(modes, atoms=None, n_confs=1000, rmsd=1.0):
 
     RMSD of the new conformation from :math:`R_0` can be calculated as
 
+
     .. math::
 
-      RMSD^k = \\sqrt{ {\\left( s \\sum_{i=1}^{m} r_i^k \\lambda^{-0.5}_i u_i  \\right)}^{2} / N } = \\frac{s}{ \\sqrt{N}} \\sqrt{ \\sum_{i=1}^{m} (r_i^k)^2 \\lambda^{-1}_i  }
+      RMSD^k = \\sqrt{  \\left[ s \\sum_{i=1}^{m} r_i^k \\lambda^{-0.5}_i u_i \\right] ^{2} / N } = \\frac{s}{ \\sqrt{N}} \\sqrt{ \\sum_{i=1}^{m} (r_i^k)^2 \\lambda^{-1}_i  }
 
 
     Average :math:`RMSD` of the generated conformations from the initial conformation is:
diff --git a/prody/proteins/fixer.py b/prody/proteins/fixer.py
index ad30aa731..ccd49aef3 100644
--- a/prody/proteins/fixer.py
+++ b/prody/proteins/fixer.py
@@ -28,7 +28,7 @@ def addMissingAtoms(infile, method='openbabel', pH=7.0, outfile=None, **kwargs):
     or PDBFixer with OpenMM. 
     
     There are also options whether to *model_residues* (default False), *remove_heterogens* 
-    (default False), *keep_waters* (default True), *overwrite* (default False).
+    (default False), *keep_waters* (default True), *overwrite* (default False), *keep_ids* (default True).
     
     :arg infile: PDB file name
     :type infile: str
@@ -41,9 +41,17 @@ def addMissingAtoms(infile, method='openbabel', pH=7.0, outfile=None, **kwargs):
             default is 'openbabel'
     :type method: str
     
-    :arg pH: pH value applyed only for PDBfixer.
+    :arg pH: pH value applied only for PDBfixer.
     :type pH: int, float
     
+    :arg model_residues: add all missing atoms from residues, applied only for PDBfixer.
+                    default is False
+    :type model_residues: bool
+    
+    :arg keep_ids: keep the original residue number, applied only for PDBfixer.
+                    default is True
+    :type keep_ids: bool
+
     Instalation of Openbabel:
     conda install -c conda-forge openbabel    
 
@@ -58,6 +66,7 @@ def addMissingAtoms(infile, method='openbabel', pH=7.0, outfile=None, **kwargs):
     remove_heterogens = kwargs.get("remove_heterogens", False)
     keep_water = kwargs.get("keep_water", True)
     overwrite = kwargs.get("overwrite", False)
+    keep_ids = kwargs.get("keep_ids", True)
 
     import os
 
@@ -70,6 +79,9 @@ def addMissingAtoms(infile, method='openbabel', pH=7.0, outfile=None, **kwargs):
     if not isinstance(keep_water, bool):
         raise TypeError('keep_water should be True or False')
 
+    if not isinstance(keep_ids, bool):
+        raise TypeError('keep_ids should be True or False')
+
     if not isinstance(overwrite, bool):
         raise TypeError('overwrite should be True or False')
     
@@ -136,7 +148,7 @@ def addMissingAtoms(infile, method='openbabel', pH=7.0, outfile=None, **kwargs):
             fixer.findMissingAtoms()
             fixer.addMissingAtoms()
             fixer.addMissingHydrogens(pH)
-            PDBFile.writeFile(fixer.topology, fixer.positions, open(outfile, 'w'))
+            PDBFile.writeFile(fixer.topology, fixer.positions, open(outfile, 'w'), keepIds=keep_ids)
             LOGGER.info("Hydrogens were added to the structure. New structure is saved as {0}.".format(outfile))
 
         except ImportError:
@@ -165,8 +177,16 @@ def fixStructuresMissingAtoms(infiles, method='openbabel', pH=7.0, outfiles=None
             'pdbfixer': PDBFixer and OpenMM
             default is 'openbabel'
     :type method: str
+
+    :arg model_residues: add all missing atoms from residues, applied only for PDBfixer.
+                    default is False
+    :type model_residues: bool
+    
+    :arg keep_ids: keep the original residue number, applied only for PDBfixer.
+                    default is True
+    :type keep_ids: bool
     
-    :arg pH: pH value applyed only for PDBfixer.
+    :arg pH: pH value applied only for PDBfixer.
     :type pH: int, float
     
     Instalation of Openbabel:
diff --git a/prody/proteins/functions.py b/prody/proteins/functions.py
index 1cf051a7c..a07dcb1c1 100644
--- a/prody/proteins/functions.py
+++ b/prody/proteins/functions.py
@@ -266,7 +266,7 @@ def showProtein(*atoms, **kwargs):
                 show = child
                 break
         if show is None:
-            show = cf.add_subplot(projection="3d")
+            show = cf.add_subplot(1,1,1,projection="3d")
         from matplotlib import colors
         cnames = dict(colors.cnames)
         wcolor = kwargs.get('water', 'red').lower()
diff --git a/prody/proteins/interactions.py b/prody/proteins/interactions.py
index 40f3dd3d6..7c49cc4c0 100644
--- a/prody/proteins/interactions.py
+++ b/prody/proteins/interactions.py
@@ -28,7 +28,7 @@
 from prody.proteins import writePDB, parsePDB
 from collections import Counter
 
-from prody.trajectory import TrajBase, Trajectory
+from prody.trajectory import TrajBase, Trajectory, Frame
 from prody.ensemble import Ensemble
 
 import multiprocessing
@@ -46,7 +46,8 @@
            'calcHydrogenBondsTrajectory', 'calcHydrophobicOverlapingAreas',
            'Interactions', 'InteractionsTrajectory', 'LigandInteractionsTrajectory',
            'calcSminaBindingAffinity', 'calcSminaPerAtomInteractions', 'calcSminaTermValues',
-           'showSminaTermValues']
+           'showSminaTermValues', 'showPairEnergy', 'checkNonstandardResidues',
+           'saveInteractionsAsDummyAtoms']
 
 
 def cleanNumbers(listContacts):
@@ -153,6 +154,160 @@ def filterInteractions(list_of_interactions, atoms, **kwargs):
     return final
 
 
+def get_energy(pair, source):
+    """Return energies based on the pairs of interacting residues (without distance criteria)
+    Taking information from tabulated_energies.txt file"""
+
+    import numpy as np
+    import importlib.resources as pkg_resources    
+    
+    aa_correction = {
+        # Histidine (His)
+        'HSD': 'HIS',   # NAMD, protonated at ND1 (HID in AMBER)
+        'HSE': 'HIS',   # NAMD, protonated at NE2 (HIE in AMBER)
+        'HSP': 'HIS',   # NAMD, doubly protonated (HIP in AMBER)
+        'HID': 'HIS',   # AMBER name, protonated at ND1
+        'HIE': 'HIS',   # AMBER name, protonated at NE2
+        'HIP': 'HIS',   # AMBER name, doubly protonated
+        'HISD': 'HIS',  # GROMACS: protonated at ND1
+        'HISE': 'HIS',  # GROMACS: protonated at NE2
+        'HISP': 'HIS',  # GROMACS: doubly protonated
+
+        # Cysteine (Cys)
+        'CYX': 'CYS',   # Cystine (disulfide bridge)
+        'CYM': 'CYS',   # Deprotonated cysteine, anion
+
+        # Aspartic acid (Asp)
+        'ASH': 'ASP',   # Protonated Asp
+        'ASPP': 'ASP',
+
+        # Glutamic acid (Glu)
+        'GLH': 'GLU',   # Protonated Glu
+        'GLUP': 'GLU',  # Protonated Glu
+
+        # Lysine (Lys)
+        'LYN': 'LYS',   # Deprotonated lysine (neutral)
+
+        # Arginine (Arg)
+        'ARN': 'ARG',   # Deprotonated arginine (rare, GROMACS)
+
+        # Tyrosine (Tyr)
+        'TYM': 'TYR',   # Deprotonated tyrosine (GROMACS)
+
+        # Serine (Ser)
+        'SEP': 'SER',   # Phosphorylated serine (GROMACS/AMBER)
+
+        # Threonine (Thr)
+        'TPO': 'THR',   # Phosphorylated threonine (GROMACS/AMBER)
+
+        # Tyrosine (Tyr)
+        'PTR': 'TYR',   # Phosphorylated tyrosine (GROMACS/AMBER)
+
+        # Non-standard names for aspartic and glutamic acids in low pH environments
+        'ASH': 'ASP',   # Protonated Asp
+        'GLH': 'GLU',   # Protonated Glu
+    }
+    
+    pair = [aa_correction.get(aa, aa) for aa in pair]    
+    
+    if PY3K:
+         file_path = pkg_resources.path('prody.proteins', 'tabulated_energies.txt')
+    else:
+        file_path = pkg_resources.resource_filename('prody.proteins', 'tabulated_energies.txt')
+
+    with open(file_path) as f:
+        data = np.loadtxt(f, dtype=str)
+    
+    sources = ["IB_nosolv", "IB_solv", "CS"]
+    aa_pairs = []
+    
+    for row in data:
+        aa_pairs.append(row[0]+row[1])
+    
+    lookup = pair[0]+pair[1]
+    
+    try:
+        data_results = data[np.where(np.array(aa_pairs)==lookup)[0]][0][2:][np.where(np.array(sources)==source)][0]
+    except TypeError:
+        raise TypeError('Please replace non-standard names of residues with standard names.')
+
+    return data_results
+
+
+def checkNonstandardResidues(atoms):
+    """Check whether the atomic structure contains non-standard residues and inform to replace the name
+    to the standard one so that non-standard residues are treated in a correct way while computing
+    interactions.
+    
+    :arg atoms: an Atomic object from which residues are selected
+    :type atoms: :class:`.Atomic`
+    """
+
+    try:
+        coords = (atoms._getCoords() if hasattr(atoms, '_getCoords') else
+                    atoms.getCoords())
+    except AttributeError:
+        try:
+            checkCoords(coords)
+        except TypeError:
+            raise TypeError('coords must be an object '
+                            'with `getCoords` method')
+    
+    amino_acids = ["ALA", "ARG", "ASN", "ASP", "CYS", "GLU", "GLN", "GLY", "HIS", "ILE", 
+                   "LEU", "LYS", "MET", "PHE", "PRO", "SER", "THR", "TRP", "TYR", "VAL"]
+
+    aa_list = atoms.select('name CA').getResnames()
+    aa_list_nr = atoms.select('name CA').getResnums()
+    nonstandard = []
+    
+    for nr_i,i in enumerate(aa_list):
+        if i not in amino_acids:
+            nonstandard.append(aa_list[nr_i] + str(aa_list_nr[nr_i]))
+    
+    if len(nonstandard) > 0:
+        LOGGER.info('There are several non-standard residues in the structure.')
+        LOGGER.info('Replace the non-standard name in the PDB file with the equivalent name from the standard one if you want to include them in the interactions.')
+        LOGGER.info("Residues: {0}.".format(' '.join(nonstandard)))
+        return True
+
+    return False
+
+
+def showPairEnergy(data, **kwargs):
+    """Return energies when a list of interactions is given. Energies will be added to each pair of residues 
+    at the last position in the list. Energy is based on the residue types and not on the distances.
+    The unit of energy is kcal/mol. The energies defined as 'IB_nosolv' (non-solvent-mediated), 'IB_solv' (solvent-mediated) 
+    are taken from [OK98]_ and 'CS' from InSty paper (under preparation). 
+    Protonation of residues is not distinguished. The protonation of residues is not distinguished. 
+    Known residues such as HSD, HSE, HIE, and HID (used in MD simulations) are treated as HIS.
+    
+    :arg data: list with interactions from calcHydrogenBonds() or other types
+    :type data: list
+    
+    :arg energy_list_type: name of the list with energies 
+                            default is 'IB_solv'
+    :type energy_list_type: 'IB_nosolv', 'IB_solv', 'CS'
+    
+    
+    .. [OK98] Keskin O., Bahar I., Badretdinov A.Y., Ptitsyn O.B., Jernigan R.L., 
+    Empirical solvet-mediated potentials hold for both intra-molecular and 
+    inter-molecular inter-residues interactions,
+    *Protein Science* **1998** 7: 2578–2586.
+    """
+    
+    if not isinstance(data, list):
+        raise TypeError('list_of_interactions must be a list of interactions.')
+
+    energy_list_type = kwargs.pop('energy_list_type', 'IB_solv')
+    
+    for i in data:
+        energy = get_energy([i[0][:3], i[3][:3]], energy_list_type)
+        i.append(float(energy))
+        
+    return data
+
+
+
 def calcHydrophobicOverlapingAreas(atoms, **kwargs):
     """Provide information about hydrophobic contacts between pairs of residues based on 
     the regsurf program. To use this function compiled hpb.so is needed.
@@ -163,7 +318,7 @@ def calcHydrophobicOverlapingAreas(atoms, **kwargs):
     :arg selection: selection string of hydrophobic residues
     :type selection: str
     
-    :arg hpb_cutoff: cutoff for hydrophobic overlaping area values
+    :arg hpb_cutoff: cutoff for hydrophobic overlapping area values
         default is 0.0
     :type hpb_cutoff: float, int
     
@@ -200,7 +355,7 @@ def calcHydrophobicOverlapingAreas(atoms, **kwargs):
         lA = [ [x[i] + str(y[i]), z[i] +'_'+ str(w[i]), ch[i]] for i in range(len(x))]
             
         output = hpb.hpb((lB,lA))
-        LOGGER.info("Hydrophobic Overlaping Areas are computed.")
+        LOGGER.info("Hydrophobic Overlapping Areas are computed.")
         output_final = [i for i in output if i[-1] >= hpb_cutoff]
         
         if cumulative_values == None:            
@@ -1000,11 +1155,11 @@ def calcHydrophobic(atoms, **kwargs):
                         non_standard works too
     :type non_standard_Hph: dict
 
-    :arg zerosHPh: zero values of hydrophobic overlaping areas included
+    :arg zerosHPh: zero values of hydrophobic overlapping areas included
         default is False
     :type zerosHPh: bool
 
-    Last value in the output corresponds to the total hydrophobic overlaping area for two residues
+    Last value in the output corresponds to the total hydrophobic overlapping area for two residues
     not only for the atoms that are included in the list. Atoms that which are listed are the closest
     between two residues and they will be inluded to draw the line in VMD_.
     
@@ -1065,7 +1220,7 @@ def calcHydrophobic(atoms, **kwargs):
     aromatic_nr = list(set(zip(atoms.aromatic.getResnums(),atoms.aromatic.getChids())))   
     aromatic = list(set(atoms.aromatic.getResnames()))
     
-    # Computing hydrophobic overlaping areas for pairs of residues:
+    # Computing hydrophobic overlapping areas for pairs of residues:
     try:
         hpb_overlaping_results = calcHydrophobicOverlapingAreas(atoms_hydrophobic, cumulative_values='pairs')
     except: 
@@ -1816,11 +1971,16 @@ def showInteractionsGraph(statistics, **kwargs):
               'HIS': 'H', 'HSD': 'H','HSE': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W', 
               'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M'}
     
-    if len(statistics[0]) != 4:
-        raise TypeError('data must be a list obtained from calcStatisticsInteractions')
+
+    if isinstance(statistics, int) or isinstance(statistics, str) or isinstance(statistics, Atomic):
+        raise TypeError('input data must be a list, use calcStatisticsInteractions to obtain statistics for a particular interaction type')
+
+    if isinstance(statistics, InteractionsTrajectory) or isinstance(statistics, Interactions):
+        raise TypeError('use calcStatisticsInteractions to obtain statistics for a particular interaction type')
+        
     else:
-        if isinstance(statistics, int) or isinstance(statistics, str):
-            raise TypeError('node_size must be a list')
+        if len(statistics[0]) != 5:
+            raise TypeError('input data must be a list obtained from calcStatisticsInteractions')
 
     code = kwargs.pop('code', None)
     if code is None:
@@ -1887,8 +2047,15 @@ def showInteractionsGraph(statistics, **kwargs):
 def calcStatisticsInteractions(data, **kwargs):
     """Return the statistics of interactions from PDB Ensemble or trajectory including:
     (1) the weight for each residue pair: corresponds to the number of counts divided by the 
-    number of frames (values >1 are obtained when residue pair creates multiple contacts); 
-    (2) average distance of interactions for each pair [in Ang] and (3) standard deviation [Ang.].
+    number of frames (values >1 are obtained when the residue pair creates multiple contacts); 
+    (2) average distance of interactions for each pair [in Ang], 
+    (3) standard deviation [Ang.],
+    (4) Energy [in kcal/mol] that is not distance dependent. Energy by default is solvent-mediated
+    from [OK98]_ ('IB_solv'). To use non-solvent-mediated entries ('IB_nosolv') from [OK98]_ or
+    solvent-mediated values obtained for InSty paper ('CS', under preparation) change 
+    `energy_list_type` parameter. 
+    If energy information is not available, please check whether the pair of residues is listed in 
+    the "tabulated_energies.txt" file, which is localized in the ProDy directory.
         
     :arg data: list with interactions from calcHydrogenBondsTrajectory() or other types
     :type data: list
@@ -1898,6 +2065,11 @@ def calcStatisticsInteractions(data, **kwargs):
         default value is 0.2 (in 20% of conformations contact appeared)
     :type weight_cutoff: int, float 
     
+    :arg energy_list_type: name of the list with energies 
+                            default is 'IB_solv'
+    :type energy_list_type: 'IB_nosolv', 'IB_solv', 'CS'
+
+    
     Example of usage: 
     >>> atoms = parsePDB('PDBfile.pdb')
     >>> dcd = Trajectory('DCDfile.dcd')
@@ -1912,6 +2084,7 @@ def calcStatisticsInteractions(data, **kwargs):
     
     interactions_list = [ (jj[0]+jj[2]+'-'+jj[3]+jj[5], jj[6]) for ii in data for jj in ii]
     weight_cutoff = kwargs.pop('weight_cutoff', 0.2)
+    energy_list_type = kwargs.pop('energy_list_type', 'IB_solv')
     
     import numpy as np
     elements = [t[0] for t in interactions_list]
@@ -1920,11 +2093,21 @@ def calcStatisticsInteractions(data, **kwargs):
     for element in elements:
         if element not in stats:
             values = [t[1] for t in interactions_list if t[0] == element]
-            stats[element] = {
-                "stddev": np.round(np.std(values),6),
-                "mean": np.round(np.mean(values),6),
-                "weight": np.round(float(len(values))/len(data), 6)
-            }
+            
+            try:
+                stats[element] = {
+                    "stddev": np.round(np.std(values),6),
+                    "mean": np.round(np.mean(values),6),
+                    "weight": np.round(float(len(values))/len(data), 6),
+                    "energy": get_energy([element.split('-')[0][:3], element.split('-')[1][:3]], energy_list_type)
+                }
+            except:
+                LOGGER.warn('energy information is not available for ', element.split('-')[0][:3], element.split('-')[1][:3])
+                stats[element] = {
+                    "stddev": np.round(np.std(values),6),
+                    "mean": np.round(np.mean(values),6),
+                    "weight": np.round(float(len(values))/len(data), 6)
+                }
 
     statistic = []
     for key, value in stats.items():
@@ -1933,7 +2116,11 @@ def calcStatisticsInteractions(data, **kwargs):
             LOGGER.info("  Average [Ang.]: {}".format(value['mean']))
             LOGGER.info("  Standard deviation [Ang.]: {0}".format(value['stddev']))
             LOGGER.info("  Weight: {0}".format(value['weight']))
-            statistic.append([key, value['weight'], value['mean'], value['stddev']])
+            try:
+                LOGGER.info("  Energy [kcal/mol]: {0}".format(value['energy']))
+                statistic.append([key, value['weight'], value['mean'], value['stddev'], value['energy']])
+            except:
+                statistic.append([key, value['weight'], value['mean'], value['stddev']])
         else: pass
     
     statistic.sort(key=lambda x: x[1], reverse=True)
@@ -2017,6 +2204,80 @@ def calcDistribution(interactions, residue1, residue2=None, **kwargs):
                 LOGGER.info(i)
 
 
+def saveInteractionsAsDummyAtoms(atoms, interactions, filename, **kwargs):
+    '''Creates a PDB file which will contain protein structure and dummy atoms that will be placed between pairs
+    of interacting residues.
+    
+    :arg atoms: an Atomic object from which residues are selected
+    :type atoms: :class:`.Atomic`
+    
+    :arg interactions: list of interactions
+    :type interactions: list
+
+    :arg filename: name of the PDB file which will contain dummy atoms and protein structure
+    :type filename: str 
+    
+    :arg RESNAME_dummy: resname of the dummy atom, use 3-letter name
+                        be default is 'DUM'
+    :type RESNAME_dummy: str '''
+
+
+    try:
+        coords = (atoms._getCoords() if hasattr(atoms, '_getCoords') else
+                    atoms.getCoords())
+    except AttributeError:
+        try:
+            checkCoords(coords)
+        except TypeError:
+            raise TypeError('coords must be an object '
+                            'with `getCoords` method')
+                            
+    RESNAME_dummy = kwargs.pop('RESNAME_dummy', 'DUM')
+    
+    def calcDUMposition(coord1, coord2):
+        midpoint = [
+            (coord1[0] + coord2[0]) / 2,
+            (coord1[1] + coord2[1]) / 2,
+            (coord1[2] + coord2[2]) / 2
+        ]
+        return midpoint
+
+    all_DUMs = []
+    atoms_ = atoms.copy()
+
+    for i in interactions:
+        if len(i[1].split('_')) <= 3:
+            res1_name = 'chain '+i[2]+' and resname '+i[0][:3]+' and resid '+i[0][3:]+' and index '+' '.join(i[1].split('_')[1:])
+            res1_coords = calcCenter(atoms.select(res1_name))  
+        
+        if len(i[1].split('_')) > 3:
+            res1_name = 'chain '+i[2]+' and resname '+i[0][:3]+' and resid '+i[0][3:]+' and index '+' '.join(i[1].split('_'))
+            res1_coords = calcCenter(atoms.select(res1_name))
+
+        if len(i[4].split('_')) <= 3:
+            res2_name = 'chain '+i[5]+' and resname '+i[3][:3]+' and resid '+i[3][3:]+' and index '+' '.join(i[4].split('_')[1:])
+            res2_coords = calcCenter(atoms.select(res2_name))
+            
+        if len(i[4].split('_')) > 3:     
+            res2_name = 'chain '+i[5]+' and resname '+i[3][:3]+' and resid '+i[3][3:]+' and index '+' '.join(i[4].split('_'))
+            res2_coords = calcCenter(atoms.select(res2_name))
+
+        all_DUMs.append(calcDUMposition(res1_coords, res2_coords))
+    
+    if all_DUMs == []:
+        LOGGER.info('Lack of interactions')
+    else:
+        LOGGER.info('Creating file with dummy atoms')
+        dummyAtoms = AtomGroup()
+        coords = array([all_DUMs], dtype=float)
+        dummyAtoms.setCoords(coords)
+        dummyAtoms.setNames([RESNAME_dummy]*len(dummyAtoms))
+        dummyAtoms.setResnums(range(1, len(dummyAtoms)+1))
+        dummyAtoms.setResnames([RESNAME_dummy]*len(dummyAtoms))
+
+        writePDB(filename, atoms_+dummyAtoms)
+
+
 def listLigandInteractions(PLIP_output, **kwargs):
     """Create a list of interactions from PLIP output created using calcLigandInteractions().
     Results can be displayed in VMD. 
@@ -2688,6 +2949,7 @@ def __init__(self, title='Unknown'):
         self._atoms = None
         self._interactions = None
         self._interactions_matrix = None
+        self._interactions_matrix_en = None
         self._hbs = None
         self._sbs = None
         self._rib = None
@@ -2767,14 +3029,33 @@ def getInteractions(self, **kwargs):
         :arg selection2: selection string
         :type selection2: str 
         
+        :arg replace: Used with selection criteria to set the new one
+                      If set to **True** the selection will be replaced by the new one.
+                      Default is **False**
+        :type replace: bool
+
         Selection:
         If we want to select interactions for the particular residue or group of residues: 
             selection='chain A and resid 1 to 50'
         If we want to study chain-chain interactions:
             selection='chain A', selection2='chain B'  """
-
+        
+        replace = kwargs.pop('replace', False)
+                
         if len(kwargs) != 0:
             results = [filterInteractions(j, self._atoms, **kwargs) for j in self._interactions]
+
+            if replace == True:
+                LOGGER.info('New interactions are set')
+                self._interactions = results           
+                self._hbs = results[0]
+                self._sbs = results[1]
+                self._rib = results[2]
+                self._piStack = results[3] 
+                self._piCat = results[4]
+                self._hps = results[5]
+                self._dibs = results[6]
+
         else: 
             results = self._interactions
         
@@ -3087,6 +3368,42 @@ def buildInteractionMatrix(self, **kwargs):
         return InteractionsMap
 
 
+    def buildInteractionMatrixEnergy(self, **kwargs):
+        """Build matrix with interaction energy comming from energy of pairs of specific residues.
+        
+        :arg energy_list_type: name of the list with energies 
+                            default is 'IB_solv'
+        :type energy_list_type: 'IB_nosolv', 'IB_solv', 'CS'
+        """
+        
+        import numpy as np
+        import matplotlib
+        import matplotlib.pyplot as plt
+        from prody.dynamics.plotting import pplot
+        
+        atoms = self._atoms   
+        interactions = self._interactions
+        energy_list_type = kwargs.pop('energy_list_type', 'IB_solv')
+
+        LOGGER.info('Calculating interactions')
+        InteractionsMap = np.zeros([atoms.select('name CA').numAtoms(),atoms.select('name CA').numAtoms()])
+        resIDs = list(atoms.select('name CA').getResnums())
+        resChIDs = list(atoms.select('name CA').getChids())
+        resIDs_with_resChIDs = list(zip(resIDs, resChIDs))
+            
+        for nr_i,i in enumerate(interactions):
+            if i != []:
+                for ii in i: 
+                    m1 = resIDs_with_resChIDs.index((int(ii[0][3:]),ii[2]))
+                    m2 = resIDs_with_resChIDs.index((int(ii[3][3:]),ii[5]))
+                    scoring = get_energy([ii[0][:3], ii[3][:3]], energy_list_type)
+                    InteractionsMap[m1][m2] = InteractionsMap[m2][m1] = InteractionsMap[m1][m2] + float(scoring) 
+
+        self._interactions_matrix_en = InteractionsMap
+        
+        return InteractionsMap
+
+
     def showInteractors(self, **kwargs):
         """Display protein residues and their number of potential interactions
         with other residues from protein structure. """
@@ -3127,23 +3444,40 @@ def saveInteractionsPDB(self, **kwargs):
         
         :arg filename: name of the PDB file which will be saved for visualization,
                      it will contain the results in occupancy column.
-        :type filename: str  """
+        :type filename: str  
+        
+        :arg energy: sum of the energy between residues
+                    default is False
+        :type energy: bool 
+        """
+        
+        energy = kwargs.pop('energy', False)
         
         if not hasattr(self, '_interactions_matrix') or self._interactions_matrix is None:
             raise ValueError('Please calculate interactions matrix first.')
 
+        if not isinstance(energy, bool):
+            raise TypeError('energy should be True or False')
+        
         import numpy as np
+        from collections import Counter
+        
+        atoms = self._atoms 
         interaction_matrix = self._interactions_matrix
-        atoms = self._atoms     
-        freq_contacts_residues = np.sum(interaction_matrix, axis=0)
+        interaction_matrix_en = self._interactions_matrix_en
         
-        from collections import Counter
-        lista_ext = []
         atoms = atoms.select("protein and noh")
+        lista_ext = []
         aa_counter = Counter(atoms.getResindices())
         calphas = atoms.select('name CA')
+        
         for i in range(calphas.numAtoms()):
-            lista_ext.extend(list(aa_counter.values())[i]*[round(freq_contacts_residues[i], 8)])
+            if energy == True:
+                matrix_en_sum = np.sum(interaction_matrix_en, axis=0)
+                lista_ext.extend(list(aa_counter.values())[i]*[round(matrix_en_sum[i], 8)]) 
+            else:
+                freq_contacts_residues = np.sum(interaction_matrix, axis=0)            
+                lista_ext.extend(list(aa_counter.values())[i]*[round(freq_contacts_residues[i], 8)])        
 
         kw = {'occupancy': lista_ext}
         if 'filename' in kwargs:
@@ -3153,6 +3487,7 @@ def saveInteractionsPDB(self, **kwargs):
             writePDB('filename', atoms, **kw)
             LOGGER.info('PDB file saved.')
 
+
     def getFrequentInteractors(self, contacts_min=3):
         """Provide a list of residues with the most frequent interactions based 
         on the following interactions:
@@ -3255,18 +3590,22 @@ def showFrequentInteractors(self, cutoff=5, **kwargs):
                 y.append(all_y[nr_ii])
 
         if SETTINGS['auto_show']:
-            matplotlib.rcParams['font.size'] = '20' 
-            fig = plt.figure(num=None, figsize=(12,6), facecolor='w')
+            matplotlib.rcParams['font.size'] = '12' 
+            fig = plt.figure(num=None, figsize=(16,5), facecolor='w')
         
         y_pos = np.arange(len(y))
         show = plt.bar(y_pos, x, align='center', alpha=0.5, color='blue')
-        plt.xticks(y_pos, y, rotation=45, fontsize=20)
-        plt.ylabel('Number of interactions')
+        plt.xticks(y_pos, y, rotation=45, fontsize=16)
+        plt.ylabel('Number of interactions', fontsize=16)
         plt.tight_layout()
 
         if SETTINGS['auto_show']:
             showFigure()
-        return show        
+            
+        dict_counts = dict(zip(y, x))
+        dict_counts_sorted = dict(sorted(dict_counts.items(), key=lambda item: item[1], reverse=True))
+            
+        return dict_counts_sorted
 
 
     def showCumulativeInteractionTypes(self, **kwargs):
@@ -3292,7 +3631,12 @@ def showCumulativeInteractionTypes(self, **kwargs):
         :type HPh: int, float
 
         :arg DiBs: score per disulfide bond
-        :type DiBs: int, float """
+        :type DiBs: int, float 
+        
+        :arg energy: sum of the energy between residues
+                    default is False
+        :type energy: bool
+        """
 
         import numpy as np
         import matplotlib
@@ -3305,87 +3649,110 @@ def showCumulativeInteractionTypes(self, **kwargs):
                'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M', 'HSE': 'H', 'HSD': 'H'}
 
         atoms = self._atoms
-
+        energy = kwargs.pop('energy', False)
+        
+        if not isinstance(energy, bool):
+            raise TypeError('energy should be True or False')
+                    
         ResNumb = atoms.select('protein and name CA').getResnums()
         ResName = atoms.select('protein and name CA').getResnames()
         ResChid = atoms.select('protein and name CA').getChids()
         ResList = [ i[0]+str(i[1])+i[2] for i in list(zip([ aa_dic[i] for i in ResName ], ResNumb, ResChid)) ]
         
-        replace_matrix = kwargs.get('replace_matrix', False)
-        matrix_all = self._interactions_matrix
-
-        HBs = kwargs.get('HBs', 1)
-        SBs = kwargs.get('SBs', 1)
-        RIB = kwargs.get('RIB', 1)
-        PiStack = kwargs.get('PiStack', 1)
-        PiCat = kwargs.get('PiCat', 1)
-        HPh = kwargs.get('HPh', 1)
-        DiBs = kwargs.get('DiBs', 1)
-    
-        matrix_hbs = self.buildInteractionMatrix(HBs=HBs, SBs=0, RIB=0,PiStack=0,PiCat=0,HPh=0,DiBs=0)
-        matrix_sbs = self.buildInteractionMatrix(HBs=0, SBs=SBs, RIB=0,PiStack=0,PiCat=0,HPh=0,DiBs=0)
-        matrix_rib = self.buildInteractionMatrix(HBs=0, SBs=0, RIB=RIB,PiStack=0,PiCat=0,HPh=0,DiBs=0)
-        matrix_pistack = self.buildInteractionMatrix(HBs=0, SBs=0, RIB=0,PiStack=PiStack,PiCat=0,HPh=0,DiBs=0)
-        matrix_picat = self.buildInteractionMatrix(HBs=0, SBs=0, RIB=0,PiStack=0,PiCat=PiCat,HPh=0,DiBs=0)
-        matrix_hph = self.buildInteractionMatrix(HBs=0, SBs=0, RIB=0,PiStack=0,PiCat=0,HPh=HPh,DiBs=0)
-        matrix_dibs = self.buildInteractionMatrix(HBs=0, SBs=0, RIB=0,PiStack=0,PiCat=0,HPh=0,DiBs=DiBs)
-
-        matrix_hbs_sum = np.sum(matrix_hbs, axis=0)
-        matrix_sbs_sum = np.sum(matrix_sbs, axis=0)
-        matrix_rib_sum = np.sum(matrix_rib, axis=0)
-        matrix_pistack_sum = np.sum(matrix_pistack, axis=0)
-        matrix_picat_sum = np.sum(matrix_picat, axis=0)
-        matrix_hph_sum = np.sum(matrix_hph, axis=0)
-        matrix_dibs_sum = np.sum(matrix_dibs, axis=0)
-
-        width = 0.8
-        fig, ax = plt.subplots(num=None, figsize=(20,6), facecolor='w')
-        matplotlib.rcParams['font.size'] = '24'
-
-        sum_matrix = np.zeros(matrix_hbs_sum.shape)
-        pplot(sum_matrix, atoms=atoms.ca)
-
-        if HBs != 0:
-            ax.bar(ResList, matrix_hbs_sum, width, color = 'blue', bottom = 0, label='HBs')
-        sum_matrix += matrix_hbs_sum
-
-        if SBs != 0:
-            ax.bar(ResList, matrix_sbs_sum, width, color = 'yellow', bottom = sum_matrix, label='SBs')
-        sum_matrix += matrix_sbs_sum
-
-        if HPh != 0:
-            ax.bar(ResList, matrix_hph_sum, width, color = 'silver', bottom = sum_matrix, label='HPh')
-        sum_matrix += matrix_hph_sum
-        
-        if RIB != 0:
-            ax.bar(ResList, matrix_rib_sum, width, color = 'red', bottom = sum_matrix, label='RIB')
-        sum_matrix += matrix_rib_sum
-
-        if PiStack != 0:
-            ax.bar(ResList, matrix_pistack_sum, width, color = 'green', bottom = sum_matrix, label='PiStack')
-        sum_matrix += matrix_pistack_sum
-
-        if PiCat != 0:
-            ax.bar(ResList, matrix_picat_sum, width, color = 'orange', bottom = sum_matrix, label='PiCat')
-        sum_matrix += matrix_picat_sum
-        
-        if DiBs != 0:
-            ax.bar(ResList, matrix_dibs_sum, width, color = 'black', bottom = sum_matrix, label='DiBs')
-        sum_matrix += matrix_dibs_sum
-
-        if replace_matrix:
-            self._interactions_matrix = np.sum([matrix_hbs, matrix_sbs, matrix_rib, matrix_pistack,
-                                                matrix_picat, matrix_hph, matrix_dibs], axis=0)
+        if energy == True:
+            matrix_en = self._interactions_matrix_en
+            matrix_en_sum = np.sum(matrix_en, axis=0)
+
+            width = 0.8
+            fig, ax = plt.subplots(num=None, figsize=(20,6), facecolor='w')
+            matplotlib.rcParams['font.size'] = '24'
+
+            ax.bar(ResNumb, matrix_en_sum, width, color='blue')
+            
+            plt.xlim([ResNumb[0]-0.5, ResNumb[-1]+0.5])
+            plt.tight_layout()    
+            plt.xlabel('Residue')
+            plt.ylabel('Cumulative Energy [kcal/mol]')
+            plt.show()
+            
+            return matrix_en_sum
+        
         else:
-            self._interactions_matrix = matrix_all
+            replace_matrix = kwargs.get('replace_matrix', False)
+            matrix_all = self._interactions_matrix
+
+            HBs = kwargs.get('HBs', 1)
+            SBs = kwargs.get('SBs', 1)
+            RIB = kwargs.get('RIB', 1)
+            PiStack = kwargs.get('PiStack', 1)
+            PiCat = kwargs.get('PiCat', 1)
+            HPh = kwargs.get('HPh', 1)
+            DiBs = kwargs.get('DiBs', 1)
+        
+            matrix_hbs = self.buildInteractionMatrix(HBs=HBs, SBs=0, RIB=0,PiStack=0,PiCat=0,HPh=0,DiBs=0)
+            matrix_sbs = self.buildInteractionMatrix(HBs=0, SBs=SBs, RIB=0,PiStack=0,PiCat=0,HPh=0,DiBs=0)
+            matrix_rib = self.buildInteractionMatrix(HBs=0, SBs=0, RIB=RIB,PiStack=0,PiCat=0,HPh=0,DiBs=0)
+            matrix_pistack = self.buildInteractionMatrix(HBs=0, SBs=0, RIB=0,PiStack=PiStack,PiCat=0,HPh=0,DiBs=0)
+            matrix_picat = self.buildInteractionMatrix(HBs=0, SBs=0, RIB=0,PiStack=0,PiCat=PiCat,HPh=0,DiBs=0)
+            matrix_hph = self.buildInteractionMatrix(HBs=0, SBs=0, RIB=0,PiStack=0,PiCat=0,HPh=HPh,DiBs=0)
+            matrix_dibs = self.buildInteractionMatrix(HBs=0, SBs=0, RIB=0,PiStack=0,PiCat=0,HPh=0,DiBs=DiBs)
+
+            matrix_hbs_sum = np.sum(matrix_hbs, axis=0)
+            matrix_sbs_sum = np.sum(matrix_sbs, axis=0)
+            matrix_rib_sum = np.sum(matrix_rib, axis=0)
+            matrix_pistack_sum = np.sum(matrix_pistack, axis=0)
+            matrix_picat_sum = np.sum(matrix_picat, axis=0)
+            matrix_hph_sum = np.sum(matrix_hph, axis=0)
+            matrix_dibs_sum = np.sum(matrix_dibs, axis=0)
+
+            width = 0.8
+            fig, ax = plt.subplots(num=None, figsize=(20,6), facecolor='w')
+            matplotlib.rcParams['font.size'] = '24'
+
+            sum_matrix = np.zeros(matrix_hbs_sum.shape)
+            pplot(sum_matrix, atoms=atoms.ca)
+
+            if HBs != 0:
+                ax.bar(ResList, matrix_hbs_sum, width, color = 'blue', bottom = 0, label='HBs')
+            sum_matrix += matrix_hbs_sum
+
+            if SBs != 0:
+                ax.bar(ResList, matrix_sbs_sum, width, color = 'yellow', bottom = sum_matrix, label='SBs')
+            sum_matrix += matrix_sbs_sum
+
+            if HPh != 0:
+                ax.bar(ResList, matrix_hph_sum, width, color = 'silver', bottom = sum_matrix, label='HPh')
+            sum_matrix += matrix_hph_sum
+            
+            if RIB != 0:
+                ax.bar(ResList, matrix_rib_sum, width, color = 'red', bottom = sum_matrix, label='RIB')
+            sum_matrix += matrix_rib_sum
 
-        ax.legend(ncol=7, loc='upper center')
-        plt.ylim([0,max(sum_matrix)+3])
-        plt.tight_layout()    
-        plt.xlabel('Residue')
-        plt.ylabel('Number of counts')
-       
-        return matrix_hbs_sum, matrix_sbs_sum, matrix_rib_sum, matrix_pistack_sum, matrix_picat_sum, matrix_hph_sum, matrix_dibs_sum 
+            if PiStack != 0:
+                ax.bar(ResList, matrix_pistack_sum, width, color = 'green', bottom = sum_matrix, label='PiStack')
+            sum_matrix += matrix_pistack_sum
+
+            if PiCat != 0:
+                ax.bar(ResList, matrix_picat_sum, width, color = 'orange', bottom = sum_matrix, label='PiCat')
+            sum_matrix += matrix_picat_sum
+            
+            if DiBs != 0:
+                ax.bar(ResList, matrix_dibs_sum, width, color = 'black', bottom = sum_matrix, label='DiBs')
+            sum_matrix += matrix_dibs_sum
+
+            if replace_matrix:
+                self._interactions_matrix = np.sum([matrix_hbs, matrix_sbs, matrix_rib, matrix_pistack,
+                                                    matrix_picat, matrix_hph, matrix_dibs], axis=0)
+            else:
+                self._interactions_matrix = matrix_all
+
+            ax.legend(ncol=7, loc='upper center')
+            plt.ylim([0,max(sum_matrix)+3])
+            plt.tight_layout()    
+            plt.xlabel('Residue')
+            plt.ylabel('Number of counts')
+        
+            return matrix_hbs_sum, matrix_sbs_sum, matrix_rib_sum, matrix_pistack_sum, matrix_picat_sum, matrix_hph_sum, matrix_dibs_sum 
       
         
 class InteractionsTrajectory(object):
@@ -3581,19 +3948,60 @@ def getInteractions(self, **kwargs):
     
         :arg selection2: selection string
         :type selection2: str 
-            
+        
+        :arg replace: Used with selection criteria to set the new one
+                      If set to **True** the selection will be replaced by the new one.
+                      Default is **False**
+        :type replace: bool
+
         Selection:
         If we want to select interactions for the particular residue or group of residues: 
             selection='chain A and resid 1 to 50'
         If we want to study chain-chain interactions:
             selection='chain A', selection2='chain B'  """
-
+        
+        replace = kwargs.pop('replace', False)
+        
         if len(kwargs) != 0:
             sele_inter = []
             for i in self._interactions_traj:
                 for nr_j,j in enumerate(i):
                     sele_inter.append(filterInteractions(i[nr_j], self._atoms, **kwargs))
+            
+            if replace == True:
+                try:
+                    trajectory = self._traj
+                    numFrames = trajectory._n_csets
+                except:
+                    # If we analyze previously saved PKL file it doesn't have dcd information
+                    # We have seven type of interactions. It will give number of frames.
+                    numFrames = int(len(sele_inter)/7)
+                    
+                self._interactions_traj = sele_inter
+                self._hbs_traj = sele_inter[0:numFrames]
+                self._sbs_traj = sele_inter[numFrames:2*numFrames]
+                self._rib_traj = sele_inter[2*numFrames:3*numFrames]
+                self._piStack_traj = sele_inter[3*numFrames:4*numFrames]
+                self._piCat_traj = sele_inter[4*numFrames:5*numFrames]
+                self._hps_traj = sele_inter[5*numFrames:6*numFrames]
+                self._dibs_traj = sele_inter[6*numFrames:7*numFrames]
+                LOGGER.info('New interactions are set')
+                
+                self._interactions_nb_traj = None
+                self._interactions_matrix_traj = None
+                
+                new_interactions_nb_traj = []
+                new_interactions_nb_traj.append([ len(i) for i in self._hbs_traj ])
+                new_interactions_nb_traj.append([ len(i) for i in self._sbs_traj ])
+                new_interactions_nb_traj.append([ len(i) for i in self._rib_traj ])
+                new_interactions_nb_traj.append([ len(i) for i in self._piStack_traj ])
+                new_interactions_nb_traj.append([ len(i) for i in self._piCat_traj ])
+                new_interactions_nb_traj.append([ len(i) for i in self._hps_traj ])
+                new_interactions_nb_traj.append([ len(i) for i in self._dibs_traj ])
+                self._interactions_nb_traj = new_interactions_nb_traj
+                
             results = sele_inter
+        
         else: 
             results = self._interactions_traj
         
diff --git a/prody/proteins/localpdb.py b/prody/proteins/localpdb.py
index ed032953f..a28cd74b5 100644
--- a/prody/proteins/localpdb.py
+++ b/prody/proteins/localpdb.py
@@ -485,8 +485,8 @@ def findPDBFiles(path, case=None, **kwargs):
     for fn in iterPDBFilenames(path, sort=True, reverse=True, **kwargs):
         fn = normpath(fn)
         pdb = splitext(split(fn)[1])[0]
-        ending = splitext(splitext(split(fn)[1])[0])[1]
-        if ending == 'gz':
+        ending = splitext(split(fn)[1])[1]
+        if ending == '.gz':
             pdb = splitext(pdb)[0]
         if len(pdb) == 7 and pdb.startswith('pdb'):
             pdb = pdb[3:]
diff --git a/prody/proteins/tabulated_energies.txt b/prody/proteins/tabulated_energies.txt
new file mode 100644
index 000000000..009c8c1b5
--- /dev/null
+++ b/prody/proteins/tabulated_energies.txt
@@ -0,0 +1,400 @@
+GLY GLY -0.27	-1.77	0.0
+GLY ALA 0.08	-2.24	0.0
+GLY VAL 0.17	-2.77	0.0
+GLY ILE 0.05	-3.5	0.0
+GLY LEU 0.33	-3.79	0.0
+GLY CYS -0.26	-3.43	0.0
+GLY MET 0.31	-3.02	0.0
+GLY PHE -0.15	-3.76	0.0
+GLY TYR -0.21	-3.34	0.0
+GLY TRP -0.31	-3.77	0.0
+GLY SER 0.08	-1.68	0.0
+GLY THR -0.11	-1.81	0.0
+GLY ASP 0.11	-1.84	0.0
+GLY ASN -0.09	-1.63	0.0
+GLY GLU 0.5	-1.58	0.0
+GLY GLN -0.11	-3.02	0.0
+GLY LYS -0.17	-1.32	0.0
+GLY ARG -0.1	-2.71	0.0
+GLY HIS 0.01	-2.47	0.0
+GLY PRO 0.13	-1.01	0.0
+ALA GLY 0.08	-2.24	0.0
+ALA ALA -0.38	-3.51	0.0
+ALA VAL -0.11	-3.86	0.0
+ALA ILE -0.08	-4.45	0.0
+ALA LEU -0.09	-5.02	0.0
+ALA CYS 0.0	-3.99	0.0
+ALA MET -0.27	-4.42	0.0
+ALA PHE 0.01	-4.43	0.0
+ALA TYR 0.0	-3.96	0.0
+ALA TRP -0.38	-4.66	0.0
+ALA SER 0.08	-2.49	0.0
+ALA THR 0.14	-2.38	0.0
+ALA ASP 0.02	-2.74	0.0
+ALA ASN 0.31	-2.06	0.0
+ALA GLU 0.2	-2.69	0.0
+ALA GLN 0.08	-3.65	0.0
+ALA LYS 0.06	-1.91	0.0
+ALA ARG 0.16	-3.26	0.0
+ALA HIS 0.01	-3.29	0.0
+ALA PRO 0.19	-1.78	0.0
+VAL GLY 0.17	-2.77	0.0
+VAL ALA -0.11	-3.86	0.0
+VAL VAL -0.49	-4.86	-1.0
+VAL ILE -0.33	-5.31	-1.2
+VAL LEU -0.48	-6.03	-1.0
+VAL CYS -0.22	-4.82	-0.4
+VAL MET -0.4	-5.16	-0.6
+VAL PHE -0.27	-5.33	-1.2
+VAL TYR 0.03	-4.54	-0.8
+VAL TRP 0.2	-4.7	-1.0
+VAL SER 0.25	-2.95	-1.0
+VAL THR 0.27	-2.87	-0.7
+VAL ASP 0.11	-3.28	0.0
+VAL ASN 0.22	-2.76	-0.5
+VAL GLU 0.3	-3.2	-0.2
+VAL GLN 0.05	-4.29	-1.1
+VAL LYS 0.4	-2.19	-1.1
+VAL ARG 0.35	-3.7	0.0
+VAL HIS -0.23	-4.14	-1.1
+VAL PRO 0.16	-2.43	0.0
+ILE GLY 0.05	-3.5	0.0
+ILE ALA -0.08	-4.45	0.0
+ILE VAL -0.33	-5.31	-1.2
+ILE ILE -0.37	-5.97	-0.6
+ILE LEU -0.5	-6.67	-1.1
+ILE CYS -0.3	-5.53	-0.4
+ILE MET 0.0	-5.37	-0.7
+ILE PHE -0.44	-6.11	-0.5
+ILE TYR -0.12	-5.39	-0.5
+ILE TRP -0.26	-5.79	-1.0
+ILE SER 0.33	-3.47	-0.9
+ILE THR 0.14	-3.61	-0.7
+ILE ASP 0.54	-3.46	-0.2
+ILE ASN 0.16	-3.42	-0.5
+ILE GLU 0.13	-3.99	-0.2
+ILE GLN 0.33	-4.65	-0.6
+ILE LYS 0.32	-2.88	-1.1
+ILE ARG 0.37	-4.29	0.0
+ILE HIS -0.02	-4.55	-1.1
+ILE PRO 0.11	-3.09	0.0
+LEU GLY 0.33	-3.79	0.0
+LEU ALA -0.09	-5.02	0.0
+LEU VAL -0.48	-6.03	-1.0
+LEU ILE -0.5	-6.67	-1.1
+LEU LEU -0.42	-7.16	-0.8
+LEU CYS -0.34	-6.13	-0.4
+LEU MET -0.3	-6.24	-0.6
+LEU PHE -0.4	-6.65	-1.2
+LEU TYR 0.08	-5.67	-1.3
+LEU TRP -0.31	-6.4	-1.0
+LEU SER 0.26	-4.12	-0.9
+LEU THR 0.04	-4.28	-0.7
+LEU ASP 0.41	-4.16	-0.2
+LEU ASN 0.23	-3.94	-0.5
+LEU GLU 0.34	-4.35	-0.2
+LEU GLN 0.17	-5.36	-1.1
+LEU LYS 0.52	-3.24	-1.1
+LEU ARG 0.21	-5.01	0.0
+LEU HIS 0.25	-4.85	-1.1
+LEU PRO 0.01	-3.75	0.0
+CYS GLY -0.26	-3.43	0.0
+CYS ALA 0.0	-3.99	0.0
+CYS VAL -0.22	-4.82	-0.4
+CYS ILE -0.3	-5.53	-0.4
+CYS LEU -0.34	-6.13	-0.4
+CYS CYS -2.38	-7.23	0.0
+CYS MET -0.08	-5.07	-0.4
+CYS PHE 0.49	-4.81	-1.3
+CYS TYR 0.1	-4.71	-1.3
+CYS TRP 0.44	4.7	-0.6
+CYS SER 0.02	-3.41	0.0
+CYS THR 0.46	-2.92	-0.8
+CYS ASP 0.06	-3.56	-0.2
+CYS ASN 1.46	-1.76	-0.9
+CYS GLU 0.72	-3.02	-0.2
+CYS GLN 0.21	-4.37	-0.9
+CYS LYS -0.5	-2.28	-0.2
+CYS ARG -0.12	-4.41	-0.2
+CYS HIS -0.64	-4.79	-0.8
+CYS PRO -0.15	-2.97	0.0
+MET GLY 0.31	-3.02	0.0
+MET ALA -0.27	-4.42	0.0
+MET VAL -0.4	-5.16	-0.6
+MET ILE 0.0	-5.37	-0.7
+MET LEU -0.3	-6.24	-0.6
+MET CYS -0.08	-5.07	-0.4
+MET MET -0.74	-5.89	-0.7
+MET PHE -0.14	-5.58	-0.8
+MET TYR -0.05	-5.01	-1.2
+MET TRP -0.2	-5.49	-1.0
+MET SER 0.25	-3.33	-0.6
+MET THR 0.19	-3.33	-1.0
+MET ASP 0.09	-3.69	0.0
+MET ASN 0.32	-3.05	-0.8
+MET GLU 0.13	-3.76	-0.2
+MET GLN 0.27	-4.46	-0.5
+MET LYS 0.61	-2.36	0.0
+MET ARG 0.3	-4.13	-0.8
+MET HIS -0.17	-4.47	-0.5
+MET PRO -0.14	-3.11	0.0
+PHE GLY -0.15	-3.76	0.0
+PHE ALA 0.01	-4.43	0.0
+PHE VAL -0.27	-5.33	-1.2
+PHE ILE -0.44	-6.11	-0.5
+PHE LEU -0.4	-6.65	-1.2
+PHE CYS 0.49	-4.81	-1.3
+PHE MET -0.14	-5.58	-0.8
+PHE PHE -0.71	-6.45	-1.2
+PHE TYR -0.07	-5.33	-1.2
+PHE TRP -0.13	-5.72	-1.2
+PHE SER -0.14	-4.01	-0.6
+PHE THR -0.02	-3.85	0.0
+PHE ASP 0.55	-3.52	0.0
+PHE ASN -0.22	-3.89	-0.2
+PHE GLU 0.09	-4.1	-0.4
+PHE GLN 0.14	-4.89	-1.7
+PHE LYS 0.57	-2.7	-0.9
+PHE ARG 0.22	-4.51	-1.1
+PHE HIS 0.79	-3.82	-0.6
+PHE PRO -0.19	-3.46	0.0
+TYR GLY -0.21	-3.34	0.0
+TYR ALA 0.0	-3.96	0.0
+TYR VAL 0.03	-4.54	-0.8
+TYR ILE -0.12	-5.39	-0.5
+TYR LEU 0.08	-5.67	-1.3
+TYR CYS 0.1	-4.71	-1.3
+TYR MET -0.05	-5.01	-1.2
+TYR PHE -0.07	-5.33	-1.2
+TYR TYR 0.2	-4.58	-1.5
+TYR TRP -0.02	-5.12	-1.2
+TYR SER -0.03	-3.43	-0.2
+TYR THR 0.0	-3.33	-0.6
+TYR ASP -0.04	-3.63	-0.6
+TYR ASN 0.04	-3.14	0.0
+TYR GLU 0.08	-3.62	-0.6
+TYR GLN -0.21	-4.75	-0.7
+TYR LYS 0.14	-2.64	-0.1
+TYR ARG -0.07	-4.32	-2.1
+TYR HIS 0.34	-3.78	-1.1
+TYR PRO -0.13	-2.92	0.0
+TRP GLY -0.31	-3.77	0.0
+TRP ALA -0.38	-4.66	0.0
+TRP VAL 0.2	-4.7	-1.0
+TRP ILE -0.26	-5.79	-1.0
+TRP LEU -0.31	-6.4	-1.0
+TRP CYS 0.44	4.7	-0.6
+TRP MET -0.2	-5.49	-1.0
+TRP PHE -0.13	-5.72	-1.2
+TRP TYR -0.02	-5.12	-1.2
+TRP TRP 0.27	-5.16	-1.2
+TRP SER 0.86	-2.86	-0.9
+TRP THR 0.18	-3.48	-0.9
+TRP ASP 0.11	-3.81	-0.5
+TRP ASN 0.2	-3.31	-1.3
+TRP GLU 0.2	-3.84	-0.5
+TRP GLN 0.01	-4.87	-1.3
+TRP LYS 1.0	-3.12	-1.2
+TRP ARG 0.04	-4.54	-1.8
+TRP HIS -0.05	-4.5	-1.7
+TRP PRO -0.85	-3.96	0.0
+SER GLY 0.08	-1.68	0.0
+SER ALA 0.08	-2.49	0.0
+SER VAL 0.25	-2.95	-1.0
+SER ILE 0.33	-3.47	-0.9
+SER LEU 0.26	-4.12	-0.9
+SER CYS 0.02	-3.41	0.0
+SER MET 0.25	-3.33	-0.6
+SER PHE -0.14	-4.01	-0.6
+SER TYR -0.03	-3.43	-0.2
+SER TRP 0.86	-2.86	-0.9
+SER SER -0.14	-2.14	-0.4
+SER THR -0.27	-2.22	-0.9
+SER ASP -0.13	-2.34	-0.2
+SER ASN -0.32	-2.13	0.2
+SER GLU -0.25	-2.57	-0.2
+SER GLN 0.1	-3.06	0.2
+SER LYS -0.5	-1.91	-1.2
+SER ARG 0.08	-2.78	-0.1
+SER HIS -0.38	-3.12	-0.5
+SER PRO -0.15	-1.56	0.0
+THR GLY -0.11	-1.81	0.0
+THR ALA 0.14	-2.38	0.0
+THR VAL 0.27	-2.87	-0.7
+THR ILE 0.14	-3.61	-0.7
+THR LEU 0.04	-4.28	-0.7
+THR CYS 0.46	-2.92	-0.8
+THR MET 0.19	-3.33	-1.0
+THR PHE -0.02	-3.85	0.0
+THR TYR 0.0	-3.33	-0.6
+THR TRP 0.18	-3.48	-0.9
+THR SER -0.27	-2.22	-0.9
+THR THR -0.22	-2.12	-0.9
+THR ASP -0.26	-2.41	-0.2
+THR ASN -0.3	-2.05	0.2
+THR GLU -0.07	-2.34	-0.5
+THR GLN -0.32	-3.44	0.2
+THR LYS -0.29	-1.64	-1.2
+THR ARG 0.38	-2.43	0.4
+THR HIS -0.05	-2.73	-0.5
+THR PRO 0.11	-1.24	0.0
+ASP GLY 0.11	-1.84	0.0
+ASP ALA 0.02	-2.74	0.0
+ASP VAL 0.11	-3.28	0.0
+ASP ILE 0.54	-3.46	-0.2
+ASP LEU 0.41	-4.16	-0.2
+ASP CYS 0.06	-3.56	-0.2
+ASP MET 0.09	-3.69	0.0
+ASP PHE 0.55	-3.52	0.0
+ASP TYR -0.04	-3.63	-0.6
+ASP TRP 0.11	-3.81	-0.5
+ASP SER -0.13	-2.34	-0.2
+ASP THR -0.26	-2.41	-0.2
+ASP ASP -0.07	-2.47	0.3
+ASP ASN -0.51	-2.5	-0.2
+ASP GLU 0.17	-2.35	0.3
+ASP GLN 0.21	-3.15	-0.2
+ASP LYS -0.87	-2.47	-1.4
+ASP ARG -0.87	-3.92	-1.4
+ASP HIS 0.0	-2.93	-0.4
+ASP PRO 0.36	-1.24	0.0
+ASN GLY -0.09	-1.63	0.0
+ASN ALA 0.31	-2.06	0.0
+ASN VAL 0.22	-2.76	-0.5
+ASN ILE 0.16	-3.42	-0.5
+ASN LEU 0.23	-3.94	-0.5
+ASN CYS 1.46	-1.76	-0.9
+ASN MET 0.32	-3.05	-0.8
+ASN PHE -0.22	-3.89	-0.2
+ASN TYR 0.04	-3.14	0.0
+ASN TRP 0.2	-3.31	-1.3
+ASN SER -0.32	-2.13	0.2
+ASN THR -0.3	-2.05	0.2
+ASN ASP -0.51	-2.5	-0.2
+ASN ASN -0.4	-1.99	-0.1
+ASN GLU -0.3	-2.43	0.1
+ASN GLN -0.05	-3.0	-0.1
+ASN LYS -0.44	-1.63	0.1
+ASN ARG -0.02	-2.67	0.2
+ASN HIS -0.52	-3.05	-1.2
+ASN PRO 0.19	-1.01	0.0
+GLU GLY 0.5	-1.58	0.0
+GLU ALA 0.2	-2.69	0.0
+GLU VAL 0.3	-3.2	-0.2
+GLU ILE 0.13	-3.99	-0.2
+GLU LEU 0.34	-4.35	-0.2
+GLU CYS 0.72	-3.02	-0.2
+GLU MET 0.13	-3.76	-0.2
+GLU PHE 0.09	-4.1	-0.4
+GLU TYR 0.08	-3.62	-0.6
+GLU TRP 0.2	-3.84	-0.5
+GLU SER -0.25	-2.57	-0.2
+GLU THR -0.07	-2.34	-0.5
+GLU ASP 0.17	-2.35	0.3
+GLU ASN -0.3	-2.43	0.1
+GLU GLU -0.2	-2.85	0.3
+GLU GLN 0.03	-3.45	0.1
+GLU LYS -1.11	-2.83	-1.0
+GLU ARG -0.87	-4.05	-1.6
+GLU HIS -0.1	-3.15	-0.4
+GLU PRO 0.01	-1.7	0.0
+GLN GLY -0.11	-3.02	0.0
+GLN ALA 0.08	-3.65	0.0
+GLN VAL 0.05	-4.29	-1.1
+GLN ILE 0.33	-4.65	-0.6
+GLN LEU 0.17	-5.36	-1.1
+GLN CYS 0.21	-4.37	-0.9
+GLN MET 0.27	-4.46	-0.5
+GLN PHE 0.14	-4.89	-1.7
+GLN TYR -0.21	-4.75	-0.7
+GLN TRP 0.01	-4.87	-1.3
+GLN SER 0.1	-3.06	0.2
+GLN THR -0.32	-3.44	0.2
+GLN ASP 0.21	-3.15	-0.2
+GLN ASN -0.05	-3.0	-0.1
+GLN GLU 0.03	-3.45	0.1
+GLN GLN -0.38	-4.71	-0.1
+GLN LYS 0.0	-2.56	0.1
+GLN ARG -0.11	-4.13	0.2
+GLN HIS -0.31	-4.2	-1.2
+GLN PRO -0.11	-2.67	0.0
+LYS GLY -0.17	-1.32	0.0
+LYS ALA 0.06	-1.91	0.0
+LYS VAL 0.4	-2.19	-1.1
+LYS ILE 0.32	-2.88	-1.1
+LYS LEU 0.52	-3.24	-1.1
+LYS CYS -0.5	-2.28	-0.2
+LYS MET 0.61	-2.36	0.0
+LYS PHE 0.57	-2.7	-0.9
+LYS TYR 0.14	-2.64	-0.1
+LYS TRP 1.0	-3.12	-1.2
+LYS SER -0.5	-1.91	-1.2
+LYS THR -0.29	-1.64	-1.2
+LYS ASP -0.87	-2.47	-1.4
+LYS ASN -0.44	-1.63	0.1
+LYS GLU -1.11	-2.83	-1.0
+LYS GLN 0.0	-2.56	0.1
+LYS LYS -0.22	-1.02	0.2
+LYS ARG 0.15	-2.11	0.0
+LYS HIS -0.01	-2.14	0.0
+LYS PRO 0.3	-0.5	0.0
+ARG GLY -0.1	-2.71	0.0
+ARG ALA 0.16	-3.26	0.0
+ARG VAL 0.35	-3.7	0.0
+ARG ILE 0.37	-4.29	0.0
+ARG LEU 0.21	-5.01	0.0
+ARG CYS -0.12	-4.41	-0.2
+ARG MET 0.3	-4.13	-0.8
+ARG PHE 0.22	-4.51	-1.1
+ARG TYR -0.07	-4.32	-2.1
+ARG TRP 0.04	-4.54	-1.8
+ARG SER 0.08	-2.78	-0.1
+ARG THR 0.38	-2.43	0.4
+ARG ASP -0.87	-3.92	-1.4
+ARG ASN -0.02	-2.67	0.2
+ARG GLU -0.87	-4.05	-1.6
+ARG GLN -0.11	-4.13	0.2
+ARG LYS 0.15	-2.11	0.0
+ARG ARG -0.26	-3.98	-0.6
+ARG HIS 0.35	-3.24	-0.4
+ARG PRO -0.17	-2.43	0.0
+HIS GLY 0.01	-2.47	0.0
+HIS ALA 0.01	-3.29	0.0
+HIS VAL -0.23	-4.14	-1.1
+HIS ILE -0.02	-4.55	-1.1
+HIS LEU 0.25	-4.85	-1.1
+HIS CYS -0.64	-4.79	-0.8
+HIS MET -0.17	-4.47	-0.5
+HIS PHE 0.79	-3.82	-0.6
+HIS TYR 0.34	-3.78	-1.1
+HIS TRP -0.05	-4.5	-1.7
+HIS SER -0.38	-3.12	-0.5
+HIS THR -0.05	-2.73	-0.5
+HIS ASP 0.0	-2.93	-0.4
+HIS ASN -0.52	-3.05	-1.2
+HIS GLU -0.1	-3.15	-0.4
+HIS GLN -0.31	-4.2	-1.2
+HIS LYS -0.01	-2.14	0.0
+HIS ARG 0.35	-3.24	-0.4
+HIS HIS 0.38	-3.08	-0.5
+HIS PRO 0.33	-1.8	0.0
+PRO GLY 0.13	-1.01	0.0
+PRO ALA 0.19	-1.78	0.0
+PRO VAL 0.16	-2.43	0.0
+PRO ILE 0.11	-3.09	0.0
+PRO LEU 0.01	-3.75	0.0
+PRO CYS -0.15	-2.97	0.0
+PRO MET -0.14	-3.11	0.0
+PRO PHE -0.19	-3.46	0.0
+PRO TYR -0.13	-2.92	0.0
+PRO TRP -0.85	-3.96	0.0
+PRO SER -0.15	-1.56	0.0
+PRO THR 0.11	-1.24	0.0
+PRO ASP 0.36	-1.24	0.0
+PRO ASN 0.19	-1.01	0.0
+PRO GLU 0.01	-1.7	0.0
+PRO GLN -0.11	-2.67	0.0
+PRO LYS 0.3	-0.5	0.0
+PRO ARG -0.17	-2.43	0.0
+PRO HIS 0.33	-1.8	0.0
+PRO PRO -0.03	-0.83	0.0
diff --git a/prody/proteins/waterbridges.py b/prody/proteins/waterbridges.py
index a247e07b5..75510b612 100644
--- a/prody/proteins/waterbridges.py
+++ b/prody/proteins/waterbridges.py
@@ -15,7 +15,18 @@
 
 from itertools import combinations
 from collections import deque
-from enum import Enum, auto
+
+from prody import PY3K
+if PY3K:
+    from enum import Enum, auto
+else:
+    from enum import Enum
+    import enum
+    from itertools import count
+    def auto(it=count()):
+        return next(it)
+    enum.auto=auto
+
 from copy import copy
 
 from prody import LOGGER, SETTINGS
@@ -263,8 +274,10 @@ def getInfoOutput(waterBridgesAtomic):
         bridgeOutput = []
 
         for atom in bridge.proteins:
-            residueInfo = f"{atom.getResname()}{atom.getResnum()}"
-            atomInfo = f"{atom.getName()}_{atom.getIndex()}"
+            residueInfo = "{0}{1}".format(atom.getResname(), 
+                                          atom.getResnum())
+            atomInfo = "{0}_{1}".format(atom.getName(),
+                                        atom.getIndex())
             chainInfo = atom.getChid()
             bridgeOutput += [residueInfo, atomInfo, chainInfo]
 
@@ -273,7 +286,9 @@ def getInfoOutput(waterBridgesAtomic):
 
         bridgeOutput += [len(bridge.waters)]
         bridgeOutput += [
-            list(map(lambda w: f"{w.getChid()}_{w.getIndex()}", bridge.waters))]
+            list(map(lambda w: "{0}_{1}".format(w.getChid(),
+                                                w.getIndex()), 
+                     bridge.waters))]
 
         output.append(bridgeOutput)
 
@@ -314,7 +329,7 @@ def getAtomicOutput(waterBridges, relations):
 
 
 def getElementsRegex(elements):
-    return f'[{"|".join(elements)}].*'
+    return '[{0}].*'.format("|".join(elements))
 
 
 def calcWaterBridges(atoms, **kwargs):
@@ -417,7 +432,7 @@ def calcWaterBridges(atoms, **kwargs):
 
     relations = RelationList(len(atoms))
     tooFarAtoms = atoms.select(
-        f'water and not within {distWR} of protein')
+        'water and not within {0} of protein'.format(distWR))
     if tooFarAtoms is None:
         consideredAtoms = atoms
     else:
@@ -434,9 +449,10 @@ def calcWaterBridges(atoms, **kwargs):
         relations[oxygen].hydrogens.append(hydrogen)
 
     proteinHydrophilic = consideredAtoms.select(
-        f'protein and name "{getElementsRegex(set(donors+acceptors))}" and within {distWR} of water')
+        'protein and name "{0}" and within {1} of water'.format(
+            getElementsRegex(set(donors+acceptors)), distWR))
 
-    proteinHydrogens = consideredAtoms.select(f'protein and hydrogen') or []
+    proteinHydrogens = consideredAtoms.select('protein and hydrogen') or []
     proteinHydroPairs = findNeighbors(
         proteinHydrophilic, DIST_COVALENT_H, proteinHydrogens) if proteinHydrogens else []
     for hydrophilic in proteinHydrophilic:
@@ -474,7 +490,8 @@ def calcWaterBridges(atoms, **kwargs):
         waterBridgesWithIndices = getUniqueElements(
             waterBridgesWithIndices, getChainBridgeTuple)
 
-    log_string = f'{len(waterBridgesWithIndices)} water bridges detected using method {method}'
+    log_string = '{0} water bridges detected using method {1}'.format(
+        len(waterBridgesWithIndices), method)
     if prefix != '':
         log_string += ' for ' + prefix
     LOGGER.info(log_string)
@@ -715,7 +732,7 @@ def getResInfo(atoms):
     chids = atoms.select('protein').getChids()
 
     for i, num in enumerate(nums):
-        dict[num] = f"{names[i]}{num}{chids[i]}"
+        dict[num] = "{names[i]}{num}{chids[i]}"
 
     return dict
 
@@ -796,7 +813,8 @@ def calcWaterBridgesStatistics(frames, trajectory, **kwargs):
     interactionCount.removeDuplicateKeys(
         lambda keys, key: (key[1], key[0]) in keys)
 
-    tableHeader = f'{"RES1":<15}{"RES2":<15}{"PERC":<10}{"DIST_AVG":<10}{"DIST_STD":<10}'
+    tableHeader = '{0:<15}{1:<15}{2:<10}{3:<10}{4:<10}'.format(
+        "RES1", "RES2", "PERC", "DIST_AVG", "DIST_STD")
     LOGGER.info(tableHeader)
     info = {}
     file = open(filename, 'w') if filename else None
@@ -819,7 +837,8 @@ def calcWaterBridgesStatistics(frames, trajectory, **kwargs):
             key1, key2 = (x, y), (y, x)
             info[key1], info[key2] = pairInfo, pairInfo
 
-        tableRow = f'{resNames[x]:<15}{resNames[y]:<15}{percentage:<10.3f}{distAvg:<10.3f}{distStd:<10.3f}'
+        tableRow = '{0:<15}{1:<15}{2:<10.3f}{3:<10.3f}{4:<10.3f}'.format(
+            resNames[x], resNames[y], percentage, distAvg, distStd)
         LOGGER.info(tableRow)
         if file:
             file.write(tableRow + '\n')
@@ -885,7 +904,7 @@ def mofifyBeta(bridgeFrames, atoms):
     atoms.setBetas(0)
     for resnum, value in residueOccurances.items():
         residueAtoms = atoms.select(
-            f'resnum {resnum}')
+            'resnum {0}'.format(resnum))
         beta = value/len(bridgeFrames)
 
         residueAtoms.setBetas(beta)
@@ -1075,7 +1094,7 @@ def calcWaterBridgesDistribution(frames, res_a, res_b=None, **kwargs):
         plt.hist(result, rwidth=0.95, density=True)
         plt.xlabel('Value')
         plt.ylabel('Probability')
-        plt.title(f'Distribution: {metric}')
+        plt.title('Distribution: {0}'.format(metric))
         if SETTINGS['auto_show']:
             showFigure()
 
@@ -1108,7 +1127,8 @@ def savePDBWaterBridges(bridges, atoms, filename):
     waterOxygens = reduceTo1D(
         bridges, lambda w: w.getIndex(), lambda b: b.waters)
     waterAtoms = atoms.select(
-        f'same residue as water within 1.6 of index {" ".join(map(str, waterOxygens))}')
+        'same residue as water within 1.6 of index {0}'.format(
+            " ".join(map(str, waterOxygens))))
 
     atomsToSave = proteinAtoms.toAtomGroup() + waterAtoms.toAtomGroup()
     return writePDB(filename, atomsToSave)
@@ -1147,21 +1167,23 @@ def saveBridgesFrame(trajectory, atoms, frameIndex, frame):
 
         waterAtoms = reduceTo1D(frame, sublistSel=lambda b: b.waters)
         waterResidues = atoms.select(
-            f'same residue as water within 1.6 of index {" ".join(map(lambda a: str(a.getIndex()), waterAtoms))}')
+            'same residue as water within 1.6 of index {0}'.format(
+                " ".join(map(lambda a: str(a.getIndex()), waterAtoms))))
 
         bridgeProteinAtoms = reduceTo1D(
             frame, lambda p: p.getResnum(), lambda b: b.proteins)
         atoms.setOccupancies(0)
-        atoms.select(
-            f'resid {" ".join(map(str, bridgeProteinAtoms))}').setOccupancies(1)
+        atoms.select('resid {0}'.format(
+            " ".join(map(str, bridgeProteinAtoms)))).setOccupancies(1)
 
         atomsToSave = atoms.select(
             'protein').toAtomGroup() + waterResidues.toAtomGroup()
 
         if trajectory:
-            writePDB(f'{filename}_{frameIndex}.pdb', atomsToSave)
+            writePDB('{0}_{1}.pdb'.format(filename, frameIndex),
+                     atomsToSave)
         else:
-            writePDB(f'{filename}_{frameIndex}.pdb',
+            writePDB('{0}_{1}.pdb'.format(filename, frameIndex),
                      atomsToSave, csets=frameIndex)
 
     if max_proc == 1:
@@ -1210,7 +1232,7 @@ def saveWaterBridges(atomicBridges, filename):
     if isInfoOutput:
         info = getWaterBridgesInfoOutput(atomicBridges)
         for frameIndex, frame in enumerate(info):
-            file.write(f'FRAME {frameIndex}\n')
+            file.write('FRAME {0}\n'.format(frameIndex))
             for bridge in frame:
                 file.write(' '.join(map(str, bridge)) + '\n')
 
diff --git a/prody/tests/proteins/test_insty.py b/prody/tests/proteins/test_insty.py
index 0fd6a47ac..fefe27258 100644
--- a/prody/tests/proteins/test_insty.py
+++ b/prody/tests/proteins/test_insty.py
@@ -12,6 +12,8 @@
 from prody.proteins.interactions import calcPiCationTrajectory, calcHydrophobicTrajectory
 from prody.proteins.interactions import calcDisulfideBondsTrajectory, calcProteinInteractions
 
+import sys
+
 class TestInteractions(unittest.TestCase):
 
     def setUp(self):
@@ -55,7 +57,7 @@ def setUp(self):
             self.data_disu = calcDisulfideBondsTrajectory(self.ATOMS)
             np.save('test_2k39_disu.npy', np.array(self.data_disu, dtype=object), allow_pickle=True)
 
-    def testAllInsteractions(self):
+    def testAllInteractions(self):
         """Test for all types of interactions."""
 
         if prody.PY3K:        
@@ -132,4 +134,28 @@ def testDisulfideBonds(self):
              assert_equal(sorted([i[-1][-1] for i in data_test if i]), sorted([i[-1][-1] for i in self.DISU_INTERACTIONS if i]),
                           'failed to get correct disulfide bonds')
         
-                     
+    def testImportHpb(self):
+
+        try:
+            import prody.proteins.hpb as hpb
+            imported_hpb = True
+        except ImportError:
+            try:
+                import hpb
+                imported_hpb = True
+            except ImportError:
+                imported_hpb = False
+            
+        if sys.version_info[1] < 11:
+            self.assertTrue(imported_hpb)
+        else:
+            self.assertFalse(imported_hpb)
+
+    @classmethod
+    def tearDownClass(cls):
+        if prody.PY3K:
+            import os
+            for filename in ['test_2k39_all.npy', 'test_2k39_hbs.npy', 'test_2k39_sbs.npy', 
+                            'test_2k39_rib.npy', 'test_2k39_PiStack.npy', 'test_2k39_PiCat.npy',
+                            'test_2k39_hph.npy', 'test_2k39_disu.npy']:
+                os.remove(filename)
diff --git a/setup.py b/setup.py
index 4c8e3a9e2..f85413b04 100644
--- a/setup.py
+++ b/setup.py
@@ -7,6 +7,8 @@
 from setuptools import setup
 from setuptools import Extension
 
+import shutil
+
 if sys.version_info[:2] < (2, 7):
     sys.stderr.write('Python 2.6 and older is not supported\n')
     sys.exit()
@@ -85,7 +87,8 @@
                     'datafiles/dcd*.dcd',
                     'datafiles/xml*.xml',
                     'datafiles/msa*',
-                    'datafiles/mmcif*cif',]
+                    'datafiles/mmcif*cif',],
+    'prody.proteins': ['tabulated_energies.txt'],
 }
 
 PACKAGE_DIR = {}
@@ -94,7 +97,15 @@
     
 from glob import glob
 tntDir = join('prody', 'utilities', 'tnt')
-hpbDir = join('prody', 'proteins', 'hpbmodule')
+hpbSoDir = join('prody', 'proteins', 'hpbmodule',
+                'hpb_Python{0}.{1}'.format(sys.version_info[0],
+                                           sys.version_info[1]))
+proteinsDir = join('prody', 'proteins')
+
+try:
+    shutil.copy(hpbSoDir + "/hpb.so", proteinsDir)
+except FileNotFoundError:
+    pass
 
 EXTENSIONS = [
     Extension('prody.dynamics.rtbtools',
@@ -164,8 +175,8 @@
 setup(
     name='ProDy',
     version=__version__,
-    author='James Krieger, She Zhang, Hongchun Li, Cihan Kaya, Ahmet Bakan, and others',
-    author_email='kriegerj@pitt.edu',
+    author='James Krieger, Karolina Mikulska-Ruminska, She Zhang, Hongchun Li, Cihan Kaya, Ahmet Bakan, and others',
+    author_email='jamesmkrieger@gmail.com',
     description='A Python Package for Protein Dynamics Analysis',
     long_description=long_description,
     url='http://www.csb.pitt.edu/ProDy',