diff --git a/PG-PuReMD/src/control.c b/PG-PuReMD/src/control.c
index 718372023d7ba6ebcac3d96991532b1019dfcc28..e477ee3367a3da09840f85faab1ebe768989c93e 100644
--- a/PG-PuReMD/src/control.c
+++ b/PG-PuReMD/src/control.c
@@ -95,6 +95,7 @@ void Read_Control_File( const char * const control_file, control_params * const
     control->cm_solver_pre_comp_droptol = 0.01;
     control->cm_solver_pre_app_type = TRI_SOLVE_PA;
     control->cm_solver_pre_app_jacobi_iters = 50;
+    control->polarization_energy_enabled = TRUE;
 
     control->T_init = 0.0;
     control->T_final = 300.;
@@ -373,6 +374,11 @@ void Read_Control_File( const char * const control_file, control_params * const
                 ival = atoi( tmp[1] );
                 control->cm_solver_pre_app_jacobi_iters = ival;
             }
+            else if ( strncmp(keyword, "include_polarization_energy", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->polarization_energy_enabled = ival;
+            }
             else if ( strncmp(tmp[0], "temp_init", MAX_LINE) == 0 )
             {
                 val = atof(tmp[1]);
diff --git a/PG-PuReMD/src/cuda/cuda_nonbonded.cu b/PG-PuReMD/src/cuda/cuda_nonbonded.cu
index c619f3eb7a2f88ff89de16d2c2ae9e14c8826d33..409c3426cdae43cd451170b017f6c5559975b71a 100644
--- a/PG-PuReMD/src/cuda/cuda_nonbonded.cu
+++ b/PG-PuReMD/src/cuda/cuda_nonbonded.cu
@@ -1023,7 +1023,7 @@ void Cuda_Compute_NonBonded_Forces( reax_system *system, control_params *control
     }
 #endif
 
-    if ( update_energy == TRUE )
+    if ( update_energy == TRUE && control->polarization_energy_enabled == TRUE )
     {
         Cuda_Compute_Polarization_Energy( system, workspace, data );
     }
diff --git a/PG-PuReMD/src/nonbonded.c b/PG-PuReMD/src/nonbonded.c
index 55b3cbefbcaad6290ffe9e44b302c9941ae008a8..089af0c1d0707a05fca969023baa346b51be9283 100644
--- a/PG-PuReMD/src/nonbonded.c
+++ b/PG-PuReMD/src/nonbonded.c
@@ -269,7 +269,10 @@ void vdW_Coulomb_Energy( reax_system const * const system,
              data->ext_press[0], data->ext_press[1], data->ext_press[2] );
 #endif
 
-    Compute_Polarization_Energy( system, data );
+    if ( control->polarization_energy_enabled == TRUE )
+    {
+        Compute_Polarization_Energy( system, data );
+    }
 }
 
 
@@ -389,7 +392,10 @@ void Tabulated_vdW_Coulomb_Energy( reax_system const * const system,
         }
     }
 
-    Compute_Polarization_Energy( system, data );
+    if ( control->polarization_energy_enabled == TRUE )
+    {
+        Compute_Polarization_Energy( system, data );
+    }
 }
 
 
diff --git a/PG-PuReMD/src/reax_types.h b/PG-PuReMD/src/reax_types.h
index 661926c15fe241e785436cb74c0d80fe6f5f60a3..0464125552278fad44128c8b632af65e54ce9b95 100644
--- a/PG-PuReMD/src/reax_types.h
+++ b/PG-PuReMD/src/reax_types.h
@@ -1585,6 +1585,8 @@ struct control_params
     /* num. of iterations used to apply preconditioner via
      * Jacobi relaxation scheme (truncated Neumann series) */
     unsigned int cm_solver_pre_app_jacobi_iters;
+    /* TRUE if polarization energy calculation is enabled, FALSE otherwise */
+    unsigned int polarization_energy_enabled;
 
     /* initial temperature of simulation, in Kelvin */
     real T_init;
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index 8803e10bac2f2604f23103512814c4bc90426d09..da6c1ffd7c769fc514c875827f8703560592d417 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -219,6 +219,10 @@ int Set_Control_Parameter( const char * const keyword,
     {
         control->cm_solver_pre_app_jacobi_iters = sstrtol( values[0], __FILE__, __LINE__ );
     }
+    else if ( strncmp(keyword, "include_polarization_energy", MAX_LINE) == 0 )
+    {
+        control->polarization_energy_enabled = sstrtol( values[0], __FILE__, __LINE__ );
+    }
     else if ( strncmp(keyword, "temp_init", MAX_LINE) == 0 )
     {
         val = sstrtod( values[0], __FILE__, __LINE__ );
@@ -530,6 +534,7 @@ void Set_Control_Defaults( reax_system * const system,
     control->cm_solver_pre_comp_droptol = 0.01;
     control->cm_solver_pre_app_type = TRI_SOLVE_PA;
     control->cm_solver_pre_app_jacobi_iters = 50;
+    control->polarization_energy_enabled = TRUE;
 
     control->molec_anal = NO_ANALYSIS;
     control->freq_molec_anal = 0;
diff --git a/sPuReMD/src/nonbonded.c b/sPuReMD/src/nonbonded.c
index 807226d532052aea191e3889115d4ea3e6ec9627..798a13c1be9972ebb51a0c35341bb5d97f527da4 100644
--- a/sPuReMD/src/nonbonded.c
+++ b/sPuReMD/src/nonbonded.c
@@ -514,7 +514,10 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
             data->press[0][0], data->press[1][1], data->press[2][2] );
 #endif
 
-    Compute_Polarization_Energy( system, control, data, workspace );
+    if ( control->polarization_energy_enabled == TRUE )
+    {
+        Compute_Polarization_Energy( system, control, data, workspace );
+    }
 }
 
 
@@ -753,7 +756,10 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
     data->E_vdW += e_vdW_total;
     data->E_Ele += e_ele_total;
 
-    Compute_Polarization_Energy( system, control, data, workspace );
+    if ( control->polarization_energy_enabled == TRUE )
+    {
+        Compute_Polarization_Energy( system, control, data, workspace );
+    }
 }
 
 
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 1220988fe0f7ccc75eb6680d07b7aa48f0c0dfe5..90c515b35f032bb6a78a0cd3fd3a0b7a34d8858e 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -1070,6 +1070,8 @@ struct control_params
     /* num. of iterations used to apply preconditioner via
      * Jacobi relaxation scheme (truncated Neumann series) */
     unsigned int cm_solver_pre_app_jacobi_iters;
+    /* TRUE if polarization energy calculation is enabled, FALSE otherwise */
+    unsigned int polarization_energy_enabled;
     /**/
     int molec_anal;
     /**/