updated average

mdapy/steinhardt_bond_orientation.py

@@ -301,6 +301,7 @@ def __init__(
         use_weight=False,
         weight=None,
         voronoi=False,
+        average=False,
     ):
         self.rc = rc
         self.nnn = nnn
@@ -367,6 +368,7 @@ def __init__(
         self.use_weight = use_weight
         self.weight = weight
         self.voronoi = voronoi
+        self.average = average
 
     @ti.func
     def _factorial(self, n: int) -> ti.f64:
@@ -472,6 +474,8 @@ def _compute(
         cglist: ti.types.ndarray(),
         inverse_box: ti.types.ndarray(element_dim=1),
         weight: ti.types.ndarray(),
+        a_qnm_r: ti.types.ndarray(),
+        a_qnm_i: ti.types.ndarray(),
     ):
         MY_EPSILON = 2.220446049250313e-15
         N = pos.shape[0]
@@ -544,6 +548,42 @@ def _compute(
                     qnm_r[i, il, m] *= facn
                     qnm_i[i, il, m] *= facn
 
+        if ti.static(self.average):
+            for i in range(N):
+                for j in range(qnm_r.shape[1]):
+                    for k in range(qnm_r.shape[2]):
+                        a_qnm_r[i, j, k] = qnm_r[i, j, k]
+                        a_qnm_i[i, j, k] = qnm_i[i, j, k]
+
+            ti.loop_config(serialize=True)
+            for i in range(N):
+                if self.voronoi:
+                    K = neighbor_number[i]
+                else:
+                    if self.nnn > 0:
+                        K = self.nnn
+                    else:
+                        K = neighbor_number[i]
+
+                N_neigh = 1
+                for jj in range(K):
+                    j = verlet_list[i, jj]
+
+                    if j >= 0:  # for voronoi neighbor
+                        for il in range(nqlist):
+                            l = qlist[il]
+                            for m in range(2 * l + 1):
+                                qnm_r[i, il, m] += a_qnm_r[j, il, m]
+                                qnm_i[i, il, m] += a_qnm_i[j, il, m]
+                        N_neigh += 1
+
+                for il in range(nqlist):
+                    l = qlist[il]
+                    for m in range(2 * l + 1):
+                        qnm_r[i, il, m] /= N_neigh
+                        qnm_i[i, il, m] /= N_neigh
+
+        for i in range(N):
             for il in range(nqlist):
                 l = qlist[il]
                 qnormfac = ti.sqrt(4 * ti.math.pi / (2 * l + 1))
@@ -614,6 +654,13 @@ def compute(self):
         qmax = self.qlist.max()
         self.qnm_r = np.zeros((self.pos.shape[0], self.nqlist, 2 * qmax + 1))
         self.qnm_i = np.zeros_like(self.qnm_r)
+        if self.average:
+            a_qnm_r = np.zeros_like(self.qnm_r)
+            a_qnm_i = np.zeros_like(self.qnm_i)
+        else:
+            a_qnm_r = np.zeros((2, 2, 2))
+            a_qnm_i = np.zeros((2, 2, 2))
+
         idxcg_count = self._get_idx(self.qlist)
         cglist = np.zeros(idxcg_count)
         if self.wlflag or self.wlhatflag:
@@ -662,6 +709,8 @@ def compute(self):
                 cglist,
                 self.inverse_box,
                 self.weight,
+                a_qnm_r,
+                a_qnm_i,
             )
         else:
             self._compute(
@@ -677,6 +726,8 @@ def compute(self):
                 cglist,
                 self.inverse_box,
                 np.zeros((2, 2)),
+                a_qnm_r,
+                a_qnm_i,
             )
         if self.old_N is not None:
             self.old_qnarray = self.qnarray.copy()
@@ -780,15 +831,16 @@ def identifySolidLiquid(self, threshold=0.7, n_bond=7):
         None,
         None,
         0.0,
-        [4, 6, 8, 10],
+        [6],
         12,
-        wlflag=True,
+        wlflag=False,
         wlhatflag=False,
+        average=True,
     )
     BO.compute()
     print(f"BO time cost: {time()-start} s.")
     print(BO.qnarray[0])
     start = time()
-    BO.identifySolidLiquid()
-    print(f"SolidLiquid time cost: {time()-start} s.")
-    print(BO.solidliquid[:10])
+    # BO.identifySolidLiquid()
+    # print(f"SolidLiquid time cost: {time()-start} s.")
+    # print(BO.solidliquid[:10])

mdapy/system.py

@@ -1170,6 +1170,7 @@ def cal_steinhardt_bond_orientation(
         wlhatflag=False,
         use_voronoi=False,
         use_weight=False,
+        average=False,
         solidliquid=False,
         max_neigh=60,
         threshold=0.7,
@@ -1284,6 +1285,7 @@ def cal_steinhardt_bond_orientation(
             use_weight,
             weight,
             use_voronoi,
+            average,
         )
         SBO.compute()
         if SBO.qnarray.shape[1] > 1:
@@ -2378,12 +2380,12 @@ def cal_lindemann_parameter(self, only_global=False):
     # nlist = neigh1.query(
     #     points, {"num_neighbors": 12, "exclude_ii": True}
     # ).toNeighborList()
-    ql = freud.order.Steinhardt(l=6, weighted=True, wl=True, wl_normalize=True)
+    ql = freud.order.Steinhardt(l=6, average=True, weighted=True)
     ql.compute(
         (box, points), neigh.nlist
     )  # neigh.nlist) #{"num_neighbors": 6, "exclude_ii": True}
 
-    print(ql.ql[:10])
+    print(ql.ql[:15])
     # print(box, points)
     system = System(box=box.Lx, pos=points.astype(np.float64))
     system.wrap_pos()
@@ -2393,9 +2395,11 @@ def cal_lindemann_parameter(self, only_global=False):
     # print(neigh.nlist[neigh.nlist[:, 0] == 0])
     # system.build_neighbor(3.1, max_neigh=23)
     system.cal_steinhardt_bond_orientation(
-        use_voronoi=True, use_weight=True, wlflag=True, wlhatflag=True
+        use_voronoi=True, average=True, use_weight=True
     )  # use_voronoi=True, use_weight=True)
-    print(system.data["ql6"].to_numpy()[:10])
+    print(system.data["ql6"].to_numpy()[:15])
+    x, y = ql.ql, system.data["ql6"].to_numpy()
+    print(np.where(x - y > 1e-3))
     # print(neigh.nlist[neigh.nlist[:, 0] == 0])
     # print(system.voro_verlet_list[0])
     # # print(neigh.nlist.distances[neigh.nlist[:, 0] == 5])