diff --git a/solventmapcreator/solvationcalculation/solvationcalculator.py b/solventmapcreator/solvationcalculation/solvationcalculator.py
index 539121287ac894eea06db55522a4f6143d6979a7..a116bf16826120dec0ce910733ab60867f0af548 100644
--- a/solventmapcreator/solvationcalculation/solvationcalculator.py
+++ b/solventmapcreator/solvationcalculation/solvationcalculator.py
@@ -31,6 +31,7 @@ First term in expression. Derivation:
 
 import logging
 import numpy as np
+import sys
 import solventmapcreator.polynomialanalysis.polynomialvaluecalculator as polyvalcalc
 
 logging.basicConfig()
@@ -77,6 +78,9 @@ def calculate_binding_energy(epsilon_i, epsilon_j, temperature, theta):
     This is in kJ/mol.
     """
     association_constant = calculate_association_constant(epsilon_i, epsilon_j, temperature)
+    #when interactions are very repulsive, then association constant is close to 0.
+    if association_constant <= sys.float_info.epsilon:
+        association_constant = sys.float_info.epsilon
     fraction_numerator = -1.0 + np.sqrt(1.0 + 8*association_constant*theta)
     fraction_denominator = 4 * association_constant * theta
     return 2 * (GAS_CONSTANT * CONVERSION_FROM_J_TO_KJ) * temperature * np.log(fraction_numerator/fraction_denominator)