diff --git a/API_reference.py b/API_reference.py
index 7368677..075d94d 100644
--- a/API_reference.py
+++ b/API_reference.py
@@ -99,6 +99,7 @@ guan.print_or_write_scattering_matrix(fermi_energy, h00, h01, length=100, on_pri
 
 # calculate topological invariant     # Source code: https://py.guanjihuan.com/source-code/calculate_topological_invariant
 chern_number = guan.calculate_chern_number_for_square_lattice(hamiltonian_function, precision=100)
+chern_number = guan.calculate_chern_number_for_square_lattice_with_Wilson_loop(hamiltonian_function, precision_of_plaquettes=10, precision_of_Wilson_loop=100)
 chern_number = guan.calculate_chern_number_for_honeycomb_lattice(hamiltonian_function, a=1, precision=300)
 wilson_loop_array = guan.calculate_wilson_loop(hamiltonian_function, k_min=-pi, k_max=pi, precision=100)
 
diff --git a/PyPI/setup.cfg b/PyPI/setup.cfg
index b515354..d182c8c 100644
--- a/PyPI/setup.cfg
+++ b/PyPI/setup.cfg
@@ -1,7 +1,7 @@
 [metadata]
 # replace with your username:
 name = guan
-version = 0.0.42
+version = 0.0.43
 author = guanjihuan
 author_email = guanjihuan@163.com
 description = An open source python package
diff --git a/PyPI/src/guan/Hamiltonian_of_finite_size_systems.py b/PyPI/src/guan/Hamiltonian_of_finite_size_systems.py
index b0955db..17e820c 100644
--- a/PyPI/src/guan/Hamiltonian_of_finite_size_systems.py
+++ b/PyPI/src/guan/Hamiltonian_of_finite_size_systems.py
@@ -3,6 +3,7 @@
 # Hamiltonian of finite size systems
 
 import numpy as np
+import guan
 
 def finite_size_along_one_direction(N, on_site=0, hopping=1, period=0):
     on_site = np.array(on_site)
@@ -106,9 +107,9 @@ def hopping_along_zigzag_direction_for_graphene(N):
     return hopping
 
 def finite_size_along_two_directions_for_graphene(N1, N2, period_1=0, period_2=0):
-    on_site = finite_size_along_one_direction(4)
-    hopping_1 = hopping_along_zigzag_direction_for_graphene(1)
+    on_site = guan.finite_size_along_one_direction(4)
+    hopping_1 = guan.hopping_along_zigzag_direction_for_graphene(1)
     hopping_2 = np.zeros((4, 4), dtype=complex)
     hopping_2[3, 0] = 1
-    hamiltonian = finite_size_along_two_directions_for_square_lattice(N1, N2, on_site, hopping_1, hopping_2, period_1, period_2)
+    hamiltonian = guan.finite_size_along_two_directions_for_square_lattice(N1, N2, on_site, hopping_1, hopping_2, period_1, period_2)
     return hamiltonian
\ No newline at end of file
diff --git a/PyPI/src/guan/calculate_topological_invariant.py b/PyPI/src/guan/calculate_topological_invariant.py
index 5059faf..e272a25 100644
--- a/PyPI/src/guan/calculate_topological_invariant.py
+++ b/PyPI/src/guan/calculate_topological_invariant.py
@@ -38,6 +38,45 @@ def calculate_chern_number_for_square_lattice(hamiltonian_function, precision=10
     chern_number = chern_number/(2*pi*1j)
     return chern_number
 
+def calculate_chern_number_for_square_lattice_with_Wilson_loop(hamiltonian_function, precision_of_plaquettes=10, precision_of_Wilson_loop=100):
+    delta = 2*pi/precision_of_plaquettes
+    chern_number = 0
+    for kx in np.arange(-pi, pi, delta):
+        for ky in np.arange(-pi, pi, delta):
+            vector_array = []
+            # line_1
+            for i0 in range(precision_of_Wilson_loop+1):
+                H_delta = hamiltonian_function(kx+delta/precision_of_Wilson_loop*i0, ky) 
+                eigenvalue, eigenvector = np.linalg.eig(H_delta)
+                vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
+                vector_array.append(vector_delta)
+            # line_2
+            for i0 in range(precision_of_Wilson_loop):
+                H_delta = hamiltonian_function(kx+delta, ky+delta/precision_of_Wilson_loop*(i0+1))  
+                eigenvalue, eigenvector = np.linalg.eig(H_delta)
+                vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
+                vector_array.append(vector_delta)
+            # line_3
+            for i0 in range(precision_of_Wilson_loop):
+                H_delta = hamiltonian_function(kx+delta-delta/precision_of_Wilson_loop*(i0+1), ky+delta)  
+                eigenvalue, eigenvector = np.linalg.eig(H_delta)
+                vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
+                vector_array.append(vector_delta)
+            # line_4
+            for i0 in range(precision_of_Wilson_loop-1):
+                H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_Wilson_loop*(i0+1))  
+                eigenvalue, eigenvector = np.linalg.eig(H_delta)
+                vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
+                vector_array.append(vector_delta)
+            Wilson_loop = 1
+            for i0 in range(len(vector_array)-1):
+                Wilson_loop = Wilson_loop*np.dot(vector_array[i0].transpose().conj(), vector_array[i0+1])
+            Wilson_loop = Wilson_loop*np.dot(vector_array[len(vector_array)-1].transpose().conj(), vector_array[0])
+            arg = np.log(np.diagonal(Wilson_loop))/1j
+            chern_number = chern_number + arg
+    chern_number = chern_number/(2*pi)
+    return chern_number
+
 def calculate_chern_number_for_honeycomb_lattice(hamiltonian_function, a=1, precision=300):
     if np.array(hamiltonian_function(0, 0)).shape==():
         dim = 1