diff --git a/src/main_gpumd/add_force.cu b/src/main_gpumd/add_force.cu
index 439ba6588..cbf1a01b0 100644
--- a/src/main_gpumd/add_force.cu
+++ b/src/main_gpumd/add_force.cu
@@ -76,7 +76,7 @@ void Add_Force::parse(const char** param, int num_param, const std::vector<Group
 
   // check the number of parameters
   if (num_param != 6 && num_param != 4) {
-    PRINT_INPUT_ERROR("electron_stop should have 5 or 3 parameters.\n");
+    PRINT_INPUT_ERROR("add_force should have 5 or 3 parameters.\n");
   }
 
   // parse grouping method
@@ -101,6 +101,12 @@ void Add_Force::parse(const char** param, int num_param, const std::vector<Group
     PRINT_INPUT_ERROR("group id should < maximum number of groups in the grouping method.\n");
   }
 
+  printf(
+    "    for atoms in group %d of grouping method %d.\n", 
+    group_id_[num_calls_], 
+    grouping_method_[num_calls_]
+  );
+
   if (num_param == 6) {
     table_length_[num_calls_] = 1;
     force_table_[num_calls_].resize(table_length_[num_calls_] * 3);
@@ -113,6 +119,9 @@ void Add_Force::parse(const char** param, int num_param, const std::vector<Group
     if (!is_valid_real(param[5], &force_table_[num_calls_][2])) {
       PRINT_INPUT_ERROR("fz should be a number.\n");
     }
+    printf("    fx = %g eV/A.\n", force_table_[num_calls_][0]);
+    printf("    fy = %g eV/A.\n", force_table_[num_calls_][1]);
+    printf("    fz = %g eV/A.\n", force_table_[num_calls_][2]);
   } else {
     std::ifstream input(param[3]);
     if (!input.is_open()) {
diff --git a/src/main_gpumd/run.cu b/src/main_gpumd/run.cu
index c406d2d17..1928a023a 100644
--- a/src/main_gpumd/run.cu
+++ b/src/main_gpumd/run.cu
@@ -17,6 +17,7 @@
 Run simulation according to the inputs in the run.in file.
 ------------------------------------------------------------------------------*/
 
+#include "add_force.cuh"
 #include "cohesive.cuh"
 #include "electron_stop.cuh"
 #include "force/force.cuh"
@@ -257,6 +258,7 @@ void Run::perform_a_run()
 #endif
 
     electron_stop.compute(time_step, atom);
+    add_force.compute(step, group, atom);
 
     integrate.compute2(time_step, double(step) / number_of_steps, group, box, atom, thermo);
 
@@ -311,6 +313,7 @@ void Run::perform_a_run()
     atom.number_of_beads);
 
   electron_stop.finalize();
+  add_force.finalize();
   integrate.finalize();
   mc.finalize();
   velocity.finalize();
@@ -461,7 +464,9 @@ void Run::parse_one_keyword(std::vector<std::string>& tokens)
     integrate.parse_move(param, num_param, group);
   } else if (strcmp(param[0], "electron_stop") == 0) {
     electron_stop.parse(param, num_param, atom.number_of_atoms, number_of_types);
-  } else if (strcmp(param[0], "mc") == 0) {
+  } else if (strcmp(param[0], "add_force") == 0) {
+    add_force.parse(param, num_param, group);
+  }else if (strcmp(param[0], "mc") == 0) {
     mc.parse_mc(param, num_param, group, atom);
   } else if (strcmp(param[0], "dftd3") == 0) {
     // nothing here; will be handled elsewhere
diff --git a/src/main_gpumd/run.cuh b/src/main_gpumd/run.cuh
index c2d319863..4670f04e1 100644
--- a/src/main_gpumd/run.cuh
+++ b/src/main_gpumd/run.cuh
@@ -19,6 +19,7 @@ class Force;
 class Integrate;
 class Measure;
 
+#include "add_force.cuh"
 #include "electron_stop.cuh"
 #include "force/force.cuh"
 #include "integrate/integrate.cuh"
@@ -68,4 +69,5 @@ private:
   MC mc;
   Measure measure;
   Electron_Stop electron_stop;
+  Add_Force add_force;
 };