diff --git a/xppbe/Model/Solutions_utils.py b/xppbe/Model/Solutions_utils.py
index b8b1a20..acff4fc 100644
--- a/xppbe/Model/Solutions_utils.py
+++ b/xppbe/Model/Solutions_utils.py
@@ -209,6 +209,9 @@ def apbs_solution(self,X):
 
         if not self.apbs_created:
             from .apbs_utils.apbs_interface import APBS
+
+            self.domain_properties['concentration'] = (self.kappa/self.ang_to_m)**2*(self.eps0*self.epsilon_2*self.kb*self.T)/(2*self.qe**2*self.Na)/1000
+            
             self.apbs_obj = APBS(self.domain_properties,self.equation,self.pqr_path)
             self.apbs_created = True
 
diff --git a/xppbe/Model/apbs_utils/apbs_interface.py b/xppbe/Model/apbs_utils/apbs_interface.py
index 657697c..b45b524 100644
--- a/xppbe/Model/apbs_utils/apbs_interface.py
+++ b/xppbe/Model/apbs_utils/apbs_interface.py
@@ -18,6 +18,7 @@ def __init__(self,domain_properties,equation,pqr_path):
         elif equation == 'nonlinear':
             equation = 'npbe'
         text = text.replace('$PBE_EQUATION',equation)
+        text = text.replace('$CONCENTRATION',str(float(domain_properties['concentration'])))
         text = text.replace('$EPSILON_1',str(float(domain_properties['epsilon_1'])))
         text = text.replace('$EPSILON_2',str(float(domain_properties['epsilon_2'])))
         text = text.replace('$TEMPERATURE',str(float(domain_properties['T'])))
diff --git a/xppbe/Model/apbs_utils/config.in b/xppbe/Model/apbs_utils/config.in
index 5de191f..7bba63a 100644
--- a/xppbe/Model/apbs_utils/config.in
+++ b/xppbe/Model/apbs_utils/config.in
@@ -9,8 +9,8 @@ elec name pbe
     mol 1
     $PBE_EQUATION
     bcfl mdh
-    ion charge 1 conc 0.125 radius 1.4
-    ion charge -1 conc 0.125 radius 1.4
+    ion charge 1 conc $CONCENTRATION radius 0.0
+    ion charge -1 conc $CONCENTRATION radius 0.0
     pdie $EPSILON_1
     sdie $EPSILON_2
     srfm smol