diff --git a/qiskit_aqua/operator.py b/qiskit_aqua/operator.py
index 861fe6315f..de976ae09e 100644
--- a/qiskit_aqua/operator.py
+++ b/qiskit_aqua/operator.py
@@ -20,6 +20,7 @@
 from functools import reduce
 import logging
 import sys
+import json
 
 import numpy as np
 from scipy import sparse as scisparse
@@ -392,17 +393,9 @@ def load_from_file(file_name):
 
         Returns:
             Operator class: the loaded operator.
-
-        Note:
-            Do we need to support complex coefficient? If so, what is the format?
         """
-        with open(file_name, 'r+') as file:
-            ham_array = file.readlines()
-        ham_array = [x.strip() for x in ham_array]
-        paulis = [[float(ham_array[2 * i + 1]), label_to_pauli(ham_array[2 * i])]
-                  for i in range(len(ham_array) // 2)]
-
-        return Operator(paulis=paulis)
+        with open(file_name, 'r') as file:
+            return Operator.load_from_dict(json.load(file))
 
     def save_to_file(self, file_name):
         """
@@ -413,10 +406,7 @@ def save_to_file(self, file_name):
 
         """
         with open(file_name, 'w') as f:
-            self._check_representation("paulis")
-            for pauli in self._paulis:
-                print("{}".format(pauli[1].to_label()), file=f)
-                print("{}".format(pauli[0]), file=f)
+            json.dump(self.save_to_dict(), f)
 
     @staticmethod
     def load_from_dict(dictionary):
@@ -472,17 +462,17 @@ def save_to_dict(self):
             dict: a dictionary contains an operator with pauli representation.
         """
         self._check_representation("paulis")
-        ret_dict = {'paulis': []}
+        ret_dict = {"paulis": []}
         for pauli in self._paulis:
-            op = {'label': pauli[1].to_label()}
+            op = {"label": pauli[1].to_label()}
             if isinstance(pauli[0], complex):
-                op['coeff'] = {'real': np.real(pauli[0]),
-                               'imag': np.imag(pauli[0])
+                op["coeff"] = {"real": np.real(pauli[0]),
+                               "imag": np.imag(pauli[0])
                                }
             else:
-                op['coeff'] = {'real': pauli[0]}
+                op["coeff"] = {"real": pauli[0]}
 
-            ret_dict['paulis'].append(op)
+            ret_dict["paulis"].append(op)
 
         return ret_dict
 
diff --git a/test/test_operator.py b/test/test_operator.py
index d4cd5c3709..9b47669846 100644
--- a/test/test_operator.py
+++ b/test/test_operator.py
@@ -19,6 +19,7 @@
 import copy
 from collections import OrderedDict
 import itertools
+import os
 
 import numpy as np
 from qiskit.tools.qi.pauli import Pauli, label_to_pauli
@@ -461,5 +462,22 @@ def test_sumbit_multiple_circutis(self):
         matrix_mode = op.eval('matrix', circuit, 'local_statevector_simulator', execute_config)[0]
         non_matrix_mode = op.eval('paulis', circuit, 'local_statevector_simulator', execute_config)[0]
 
+    def test_load_from_file(self):
+        paulis = ['IXYZ', 'XXZY', 'IIZZ', 'XXYY', 'ZZXX', 'YYYY']
+        coeffs = [0.2 + -1j * 0.8, 0.6 + -1j * 0.6, 0.8 + -1j * 0.2,
+                    -0.2 + -1j * 0.8, -0.6 - -1j * 0.6, -0.8 - -1j * 0.2]
+        op = Operator(paulis=[])
+        for coeff, pauli in zip(coeffs, paulis):
+            pauli_term = [coeff, label_to_pauli(pauli)]
+            op += Operator(paulis=[pauli_term])
+
+        op.save_to_file('temp_op.json')
+        load_op = Operator.load_from_file('temp_op.json')
+
+        self.assertTrue(os.path.exists('temp_op.json'))
+        self.assertEqual(op, load_op)
+
+        os.remove('temp_op.json')
+
 if __name__ == '__main__':
     unittest.main()