diff --git a/example/2band.atom.py b/example/2band.atom.py
index 68d7cce..2f1628f 100644
--- a/example/2band.atom.py
+++ b/example/2band.atom.py
@@ -95,7 +95,8 @@
 # G^{(2)}(i\omega;i\nu,i\nu') #
 ###############################
 
-G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params )
+# G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params )
+G2_iw = ed.G2_iw_freq_box( vec_four_indices=[(('up',0), ('dn',0), ('dn',1), ('up',1)),], **common_g2_params )[0]
 print type(G2_iw)
 print G2_iw.shape
 
diff --git a/example/slater.py b/example/slater.py
index cb975b1..ff227cb 100644
--- a/example/slater.py
+++ b/example/slater.py
@@ -145,7 +145,8 @@
 # G^{(2)}(i\omega;i\nu,i\nu') #
 ###############################
 
-G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params )
+# G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params )
+G2_iw = ed.G2_iw_freq_box( vec_four_indices=[(('up',0), ('dn',0), ('dn',1), ('up',1)),], **common_g2_params )[0]
 
 if mpi.is_master_node():
     print type(G2_iw)
diff --git a/test/python/anderson_g2_matsubara.py b/test/python/anderson_g2_matsubara.py
index a8ef2e5..26c9785 100644
--- a/test/python/anderson_g2_matsubara.py
+++ b/test/python/anderson_g2_matsubara.py
@@ -72,15 +72,33 @@
                     'n_b' : g2_n_wb, }
 
 
-# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order
-G2_iw_ph_uuuu = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params )
-G2_iw_ph_dddd = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params )
-G2_iw_ph_uudd = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params )
-G2_iw_ph_dduu = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params )
-
-# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order
-G2_iw_ph_uddu = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params )
-G2_iw_ph_duud = ed.G2_iw_legacy( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params )
+# # Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order
+# G2_iw_ph_uuuu = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params )
+# G2_iw_ph_dddd = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params )
+# G2_iw_ph_uudd = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params )
+# G2_iw_ph_dduu = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params )
+#
+# # Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order
+# G2_iw_ph_uddu = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params )
+# G2_iw_ph_duud = ed.G2_iw_legacy( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params )
+
+
+four_indices = []
+four_indices.append( (('up',0), ('up',0), ('up',0), ('up',0)) )  # uuuu
+four_indices.append( (('dn',0), ('dn',0), ('dn',0), ('dn',0)) )  # dddd
+four_indices.append( (('up',0), ('up',0), ('dn',0), ('dn',0)) )  # uudd
+four_indices.append( (('dn',0), ('dn',0), ('up',0), ('up',0)) )  # dduu
+four_indices.append( (('up',0), ('dn',0), ('dn',0), ('up',0)) )  # uddu
+four_indices.append( (('dn',0), ('up',0), ('up',0), ('dn',0)) )  # duud
+
+G2_iw_ph = ed.G2_iw_freq_box( vec_four_indices=four_indices, **common_g2_params )
+G2_iw_ph_uuuu = G2_iw_ph[0]
+G2_iw_ph_dddd = G2_iw_ph[1]
+G2_iw_ph_uudd = G2_iw_ph[2]
+G2_iw_ph_dduu = G2_iw_ph[3]
+G2_iw_ph_uddu = G2_iw_ph[4]
+G2_iw_ph_duud = G2_iw_ph[5]
+
 
 # Compute G^{(2),pp}(i\omega;i\nu,i\nu'), AABB block order
 # G2_iw_inu_inup_pp_AABB = ed.G2_iw_inu_inup(channel = "PP",
diff --git a/test/python/wick.py b/test/python/wick.py
index 858ff72..d2a3c6d 100644
--- a/test/python/wick.py
+++ b/test/python/wick.py
@@ -81,19 +81,43 @@
                     'n_f' : g2_n_wf,
                     'n_b' : g2_n_wb, }
 
-G2_ph_uuuu = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params )
-G2_ph_dddd = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params )
-G2_ph_uudd = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params )
-G2_ph_dduu = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params )
-G2_ph_uddu = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params )
-G2_ph_duud = ed.G2_iw_legacy( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params )
-
-G2_ph_uuuu_wick = G2_ph_uuuu.copy()
-G2_ph_dddd_wick = G2_ph_dddd.copy()
-G2_ph_uudd_wick = G2_ph_uudd.copy()
-G2_ph_dduu_wick = G2_ph_dduu.copy()
-G2_ph_uddu_wick = G2_ph_uddu.copy()
-G2_ph_duud_wick = G2_ph_duud.copy()
+# G2_ph_uuuu = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params )
+# G2_ph_dddd = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params )
+# G2_ph_uudd = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params )
+# G2_ph_dduu = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params )
+# G2_ph_uddu = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params )
+# G2_ph_duud = ed.G2_iw_legacy( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params )
+
+four_indices = []
+four_indices.append( (('up',0), ('up',0), ('up',0), ('up',0)) )  # uuuu
+four_indices.append( (('dn',0), ('dn',0), ('dn',0), ('dn',0)) )  # dddd
+four_indices.append( (('up',0), ('up',0), ('dn',0), ('dn',0)) )  # uudd
+four_indices.append( (('dn',0), ('dn',0), ('up',0), ('up',0)) )  # dduu
+four_indices.append( (('up',0), ('dn',0), ('dn',0), ('up',0)) )  # uddu
+four_indices.append( (('dn',0), ('up',0), ('up',0), ('dn',0)) )  # duud
+
+G2_ph = ed.G2_iw_freq_box( vec_four_indices=four_indices, **common_g2_params )
+G2_ph_uuuu = G2_ph[0]
+G2_ph_dddd = G2_ph[1]
+G2_ph_uudd = G2_ph[2]
+G2_ph_dduu = G2_ph[3]
+G2_ph_uddu = G2_ph[4]
+G2_ph_duud = G2_ph[5]
+
+# G2_ph_uuuu_wick = G2_ph_uuuu.copy()
+# G2_ph_dddd_wick = G2_ph_dddd.copy()
+# G2_ph_uudd_wick = G2_ph_uudd.copy()
+# G2_ph_dduu_wick = G2_ph_dduu.copy()
+# G2_ph_uddu_wick = G2_ph_uddu.copy()
+# G2_ph_duud_wick = G2_ph_duud.copy()
+
+g2_shape = (g2_n_wb, 2*g2_n_wf, 2*g2_n_wf)
+G2_ph_uuuu_wick = np.zeros(g2_shape, dtype=np.complex)
+G2_ph_dddd_wick = np.zeros(g2_shape, dtype=np.complex)
+G2_ph_dduu_wick = np.zeros(g2_shape, dtype=np.complex)
+G2_ph_uudd_wick = np.zeros(g2_shape, dtype=np.complex)
+G2_ph_uddu_wick = np.zeros(g2_shape, dtype=np.complex)
+G2_ph_duud_wick = np.zeros(g2_shape, dtype=np.complex)
 
 G = lambda s, i: G_iw[s].data[i + n_iw, 0, 0]