diff --git a/solventmapcreator/solvationcalculation/solvationcalculator.py b/solventmapcreator/solvationcalculation/solvationcalculator.py
index d7f12d7b7adf4412f876c7b1d6ae426a23d73978..e085aa94accadff855a75d57ad8674ecbc75ed33 100644
--- a/solventmapcreator/solvationcalculation/solvationcalculator.py
+++ b/solventmapcreator/solvationcalculation/solvationcalculator.py
@@ -11,8 +11,21 @@ where R is gas constant,
 T is temperature,
 K is association constant,
 theta is fractional occupancy of the phase,
-S1, S2 are the solvation energy of SSIP 1 and 2 respectively.
+S1, S2 are the solvation energy of SSIP 1 and 2 respectively. This is just
+the binding energy contribution, \delta G_{b}.
 
+First term in expression. Derivation:
+    Consider SSIPs i and j that interact, K_{ij}!=1, and assume self interaction is zero.
+    
+    [i] = [i_{free}] + K[i_{free}][j_{free}]
+    
+    and noting that [i_{free}] = [j_{free}]
+    probability that SSIP is unbound is:
+        P_{f} = [i_{free}]/[i] = -1 + sqrt(1 + 4K*theta)]/[2K*theta]
+    so the correction per SSIP is 
+     R* T * ln( [-1 + sqrt(1 + 4K*theta)]/[2K*theta])
+    Since both SSIPs have a probility of being free then this becomes
+    2*R * T * ln( [-1 + sqrt(1 + 4K*theta)]/[2K*theta]) 
 @author: mark
 """