diff --git a/examples/nnp-predict/H2O_RPBE-D3_4G/input.nn b/examples/nnp-predict/H2O_RPBE-D3_4G/input.nn
index 8461ea8a7..4e337eac4 100644
--- a/examples/nnp-predict/H2O_RPBE-D3_4G/input.nn
+++ b/examples/nnp-predict/H2O_RPBE-D3_4G/input.nn
@@ -26,6 +26,7 @@ random_seed 12346
 
 ewald_truncation_error_method 1
 ewald_prec 1.e-5 0.36
+four_pi_epsilon 1.0                             # Value of 4*pi*epsilon (1/Coulomb constant)
 
 ################################################################################
 ### NN structure of the short-range NN
diff --git a/examples/nnp-train/H2O_RPBE-D3_4G/input.nn b/examples/nnp-train/H2O_RPBE-D3_4G/input.nn
index 00794e72c..775d81058 100644
--- a/examples/nnp-train/H2O_RPBE-D3_4G/input.nn
+++ b/examples/nnp-train/H2O_RPBE-D3_4G/input.nn
@@ -25,6 +25,7 @@ random_seed 12346
 
 ewald_truncation_error_method 1
 ewald_prec 1.e-5 0.36
+four_pi_epsilon 1.0                             # Value of 4*pi*epsilon (1/Coulomb constant)
 
 ################################################################################
 ### NN structure of the short-range NN
diff --git a/src/libnnp/Mode.cpp b/src/libnnp/Mode.cpp
index cff81e515..d12ab363c 100644
--- a/src/libnnp/Mode.cpp
+++ b/src/libnnp/Mode.cpp
@@ -47,6 +47,7 @@ Mode::Mode() : nnpType                   (NNPType::HDNNP_2G),
                convEnergy                (1.0              ),
                convLength                (1.0              ),
                convCharge                (1.0              ),
+               fourPiEps                 (1.0              ),
                ewaldSetup                {}
 {
 }
@@ -483,6 +484,12 @@ void Mode::setupElectrostatics(bool   initialHardness,
                      ewaldSetup.getMaxCharge());
     log << "\n";
 
+    // 4 * pi * epsilon
+    if (settings.keywordExists("four_pi_epsilon"))
+        fourPiEps = atof(settings["four_pi_epsilon"].c_str());
+    if (normalize)
+        fourPiEps *= pow(convCharge, 2) / (convLength * convEnergy);
+
     // Screening function.
     if (settings.keywordExists("screen_electrostatics"))
     {
@@ -1753,12 +1760,6 @@ void Mode::chargeEquilibration( Structure& structure,
     VectorXd hardness(numElements);
     MatrixXd gammaSqrt2(numElements, numElements);
     VectorXd sigmaSqrtPi(numElements);
-    // In case of natural units and no normalization
-    double fourPiEps = 1;
-    // In case of natural units but with normalization. Other units currently
-    // not supported.
-    if (normalize)
-        fourPiEps = pow(convCharge, 2) / (convLength * convEnergy);
 
     for (size_t i = 0; i < numElements; ++i)
     {
diff --git a/src/libnnp/Mode.h b/src/libnnp/Mode.h
index d4854a5d7..cb349a557 100644
--- a/src/libnnp/Mode.h
+++ b/src/libnnp/Mode.h
@@ -637,6 +637,7 @@ class Mode
     double                     convEnergy;
     double                     convLength;
     double                     convCharge;
+    double                     fourPiEps;
     EwaldSetup                 ewaldSetup;
     settings::Settings         settings;
     SymFnc::ScalingType        scalingType;
diff --git a/src/libnnp/Settings.cpp b/src/libnnp/Settings.cpp
index d63c64e93..6efee6380 100644
--- a/src/libnnp/Settings.cpp
+++ b/src/libnnp/Settings.cpp
@@ -66,6 +66,7 @@ map<string, shared_ptr<settings::Key>> const createKnownKeywordsMap()
     m["ewald_truncation_error_method"      ] = "";
     m["ewald_prec"                         ] = "";
     m["screen_electrostatics"              ] = "";
+    m["four_pi_epsilon"                    ] = "";
 
     // Training keywords.
     m["random_seed"                        ] = "";
diff --git a/src/libnnp/Structure.cpp b/src/libnnp/Structure.cpp
index 1624fa6f3..0845beb53 100644
--- a/src/libnnp/Structure.cpp
+++ b/src/libnnp/Structure.cpp
@@ -744,7 +744,7 @@ double Structure::calculateElectrostaticEnergy(
 }
 
 double Structure::calculateScreeningEnergy(
-                                        Eigen::MatrixXd          gammaSqrt2,
+                                        MatrixXd                 gammaSqrt2,
                                         VectorXd                 sigmaSqrtPi,
                                         ScreeningFunction const& fs,
                                         double const             fourPiEps)