Compute all order parameters of site n.

        Args:
            structure (Structure): input structure.
            n (int): index of site in input structure,
                for which OPs are to be
                calculated.  Note that we do not use the sites iterator
                here, but directly access sites via struct[index].
            indices_neighs ([int]): list of indices of those neighbors
                in Structure object
                structure that are to be considered for OP computation.
                This optional argument overwrites the way neighbors are
                to be determined as defined in the constructor (i.e.,
                Voronoi coordination finder via negative cutoff radius
                vs constant cutoff radius if cutoff was positive).
                We do not use information about the underlying
                structure lattice if the neighbor indices are explicitly
                provided.  This has two important consequences.  First,
                the input Structure object can, in fact, be a
                simple list of Site objects.  Second, no nearest images
                of neighbors are determined when providing an index list.
                Note furthermore that this neighbor
                determination type ignores the optional target_spec
                argument.
            tol (float): threshold of weight
                (= solid angle / maximal solid angle)
                to determine if a particular pair is
                considered neighbors; this is relevant only in the case
                when Voronoi polyhedra are used to determine coordination
            target_spec (Specie): target species to be considered
                when calculating the order
                parameters of site n; None includes all species of input
                structure.

        Returns:
            [floats]: representing order parameters.  Should it not be
            possible to compute a given OP for a conceptual reason, the
            corresponding entry is None instead of a float.  For Steinhardt
            et al.'s bond orientational OPs and the other geometric OPs
            ("tet", "oct", "bcc", etc.),
            this can happen if there is a single
            neighbor around site n in the structure because that
            does not permit calculation of angles between multiple
            neighbors.

def get_order_parameters(self, structure, n, indices_neighs=None, \
                             tol=0.0, target_spec=None):

        """
        Compute all order parameters of site n.

        Args:
            structure (Structure): input structure.
            n (int): index of site in input structure,
                for which OPs are to be
                calculated.  Note that we do not use the sites iterator
                here, but directly access sites via struct[index].
            indices_neighs ([int]): list of indices of those neighbors
                in Structure object
                structure that are to be considered for OP computation.
                This optional argument overwrites the way neighbors are
                to be determined as defined in the constructor (i.e.,
                Voronoi coordination finder via negative cutoff radius
                vs constant cutoff radius if cutoff was positive).
                We do not use information about the underlying
                structure lattice if the neighbor indices are explicitly
                provided.  This has two important consequences.  First,
                the input Structure object can, in fact, be a
                simple list of Site objects.  Second, no nearest images
                of neighbors are determined when providing an index list.
                Note furthermore that this neighbor
                determination type ignores the optional target_spec
                argument.
            tol (float): threshold of weight
                (= solid angle / maximal solid angle)
                to determine if a particular pair is
                considered neighbors; this is relevant only in the case
                when Voronoi polyhedra are used to determine coordination
            target_spec (Specie): target species to be considered
                when calculating the order
                parameters of site n; None includes all species of input
                structure.

        Returns:
            [floats]: representing order parameters.  Should it not be
            possible to compute a given OP for a conceptual reason, the
            corresponding entry is None instead of a float.  For Steinhardt
            et al.'s bond orientational OPs and the other geometric OPs
            ("tet", "oct", "bcc", etc.),
            this can happen if there is a single
            neighbor around site n in the structure because that
            does not permit calculation of angles between multiple
            neighbors.
        """

        # Do error-checking and initialization.
        if n < 0:
            raise ValueError("Site index smaller zero!")
        if n >= len(structure):
            raise ValueError("Site index beyond maximum!")
        if indices_neighs is not None:
            for index in indices_neighs:
                if index >= len(structure):
                    raise ValueError("Neighbor site index beyond maximum!")
        if tol < 0.0:
            raise ValueError("Negative tolerance for weighted solid angle!")

        left_of_unity = 1.0 - 1.0e-12
        # The following threshold has to be adapted to non-Angstrom units.
        very_small = 1.0e-12
        fac_bcc = 1.0 / exp(-0.5)

        # Find central site and its neighbors.
        # Note that we adopt the same way of accessing sites here as in
        # VoronoiNN; that is, not via the sites iterator.
        centsite = structure[n]
        if indices_neighs is not None:
            neighsites = [structure[index] for index in indices_neighs]
        elif self._voroneigh:
            vnn = VoronoiNN(tol=tol, targets=target_spec)
            neighsites = vnn.get_nn(structure, n)
        else:
            # Structure.get_sites_in_sphere --> also other periodic images
            neighsitestmp = [i[0] for i in structure.get_sites_in_sphere(
                centsite.coords, self._cutoff)]
            neighsites = []
            if centsite not in neighsitestmp:
                raise ValueError("Could not find center site!")
            else:
                neighsitestmp.remove(centsite)
            if target_spec is None:
                neighsites = list(neighsitestmp)
            else:
                neighsites[:] = [site for site in neighsitestmp \
                                 if site.specie.symbol == target_spec]
        nneigh = len(neighsites)
        self._last_nneigh = nneigh

        # Prepare angle calculations, if applicable.
        rij = []
        rjk = []
        rijnorm = []
        rjknorm = []
        dist = []
        distjk_unique = []
        distjk = []
        centvec = centsite.coords
        if self._computerijs:
            for j, neigh in enumerate(neighsites):
                rij.append((neigh.coords - centvec))
                dist.append(np.linalg.norm(rij[j]))
                rijnorm.append((rij[j] / dist[j]))
        if self._computerjks:
            for j, neigh in enumerate(neighsites):
                rjk.append([])
                rjknorm.append([])
                distjk.append([])
                kk = 0
                for k in range(len(neighsites)):
                    if j != k:
                        rjk[j].append(neighsites[k].coords - neigh.coords)
                        distjk[j].append(np.linalg.norm(rjk[j][kk]))
                        if k > j:
                            distjk_unique.append(distjk[j][kk])
                        rjknorm[j].append(rjk[j][kk] / distjk[j][kk])
                        kk = kk + 1
        # Initialize OP list and, then, calculate OPs.
        ops = [0.0 for t in self._types]
        # norms = [[[] for j in range(nneigh)] for t in self._types]

        # First, coordination number and distance-based OPs.
        for i, t in enumerate(self._types):
            if t == "cn":
                ops[i] = nneigh / self._params[i]['norm']
            elif t == "sgl_bd":
                dist_sorted = sorted(dist)
                if len(dist_sorted) == 1:
                    ops[i] = 1.0
                elif len(dist_sorted) > 1:
                    ops[i] = 1.0 - dist_sorted[0] / dist_sorted[1]

        # Then, bond orientational OPs based on spherical harmonics
        # according to Steinhardt et al., Phys. Rev. B, 28, 784-805, 1983.
        if self._boops:
            thetas = []
            phis = []
            for j, vec in enumerate(rijnorm):

                # z is North pole --> theta between vec and (0, 0, 1)^T.
                # Because vec is normalized, dot product is simply vec[2].
                thetas.append(acos(max(-1.0, min(vec[2], 1.0))))
                tmpphi = 0.0

                # Compute phi only if it is not (almost) perfectly
                # aligned with z-axis.
                if -left_of_unity < vec[2] < left_of_unity:
                    # x is prime meridian --> phi between projection of vec
                    # into x-y plane and (1, 0, 0)^T
                    tmpphi = acos(max(-1.0, min(vec[0] / (sqrt(
                        vec[0] * vec[0] + vec[1] * vec[1])), 1.0)))
                    if vec[1] < 0.0:
                        tmpphi = -tmpphi
                phis.append(tmpphi)

            # Note that None flags that we have too few neighbors
            # for calculating BOOPS.
            for i, t in enumerate(self._types):
                if t == "q2":
                    ops[i] = self.get_q2(thetas, phis) if len(
                        thetas) > 0 else None
                elif t == "q4":
                    ops[i] = self.get_q4(thetas, phis) if len(
                        thetas) > 0 else None
                elif t == "q6":
                    ops[i] = self.get_q6(thetas, phis) if len(
                        thetas) > 0 else None

        # Then, deal with the Peters-style OPs that are tailor-made
        # to recognize common structural motifs
        # (Peters, J. Chem. Phys., 131, 244103, 2009;
        #  Zimmermann et al., J. Am. Chem. Soc., under revision, 2015).
        if self._geomops:
            gaussthetak = [0.0 for t in self._types]  # not used by all OPs
            qsptheta = [[[] for j in range(nneigh)] for t in self._types]
            norms = [[[] for j in range(nneigh)] for t in self._types]
            ipi = 1.0 / pi
            piover2 = pi / 2.0
            tetangoverpi = acos(-1.0 / 3.0) * ipi  # xxx: delete
            itetangminuspihalfoverpi = 1.0 / (tetangoverpi - 0.5)
            onethird = 1.0 / 3.0
            twothird = 2.0 / 3.0
            for j in range(nneigh):  # Neighbor j is put to the North pole.
                zaxis = rijnorm[j]
                kc = 0
                for k in range(nneigh):  # From neighbor k, we construct
                    if j != k:  # the prime meridian.
                        for i in range(len(self._types)):
                            qsptheta[i][j].append(0.0)
                            norms[i][j].append(0)
                        tmp = max(
                            -1.0, min(np.inner(zaxis, rijnorm[k]), 1.0))
                        thetak = acos(tmp)
                        xaxis = gramschmidt(rijnorm[k], zaxis)
                        if np.linalg.norm(xaxis) < very_small:
                            flag_xaxis = True
                        else:
                            xaxis = xaxis / np.linalg.norm(xaxis)
                            flag_xaxis = False
                        if self._comp_azi:
                            flag_yaxis = True
                            yaxis = np.cross(zaxis, xaxis)
                            if np.linalg.norm(yaxis) > very_small:
                                yaxis = yaxis / np.linalg.norm(yaxis)
                                flag_yaxis = False

                        # Contributions of j-i-k angles, where i represents the
                        # central atom and j and k two of the neighbors.
                        for i, t in enumerate(self._types):
                            if t in ["bent", "sq_pyr_legacy"]:
                                tmp = self._params[i]['IGW_TA'] * (
                                        thetak * ipi - self._params[i]['TA'])
                                qsptheta[i][j][kc] += exp(-0.5 * tmp * tmp)
                                norms[i][j][kc] += 1
                            elif t in ["tri_plan", "tri_plan_max", "tet",
                                       "tet_max"]:
                                tmp = self._params[i]['IGW_TA'] * (
                                        thetak * ipi - self._params[i]['TA'])
                                gaussthetak[i] = exp(-0.5 * tmp * tmp)
                                if t in ["tri_plan_max", "tet_max"]:
                                    qsptheta[i][j][kc] += gaussthetak[i]
                                    norms[i][j][kc] += 1
                            elif t in ["T", "tri_pyr", "sq_pyr", "pent_pyr",
                                       "hex_pyr"]:
                                tmp = self._params[i]['IGW_EP'] * (
                                            thetak * ipi - 0.5)
                                qsptheta[i][j][kc] += exp(-0.5 * tmp * tmp)
                                norms[i][j][kc] += 1
                            elif t in ["sq_plan", "oct", "oct_legacy",
                                       "cuboct", "cuboct_max"]:
                                if thetak >= self._params[i]['min_SPP']:
                                    tmp = self._params[i]['IGW_SPP'] * (
                                            thetak * ipi - 1.0)
                                    qsptheta[i][j][kc] += (
                                            self._params[i]['w_SPP'] *
                                            exp(-0.5 * tmp * tmp))
                                    norms[i][j][kc] += self._params[i]['w_SPP']
                            elif t in ["see_saw_rect", "tri_bipyr", "sq_bipyr",
                                       "pent_bipyr", "hex_bipyr", "oct_max",
                                       "sq_plan_max", "hex_plan_max"]:
                                if thetak < self._params[i]['min_SPP']:
                                    tmp = self._params[i]['IGW_EP'] * (
                                            thetak * ipi - 0.5) if t != "hex_plan_max" else \
                                        self._params[i]['IGW_TA'] * (
                                                fabs(thetak * ipi - 0.5) -
                                                self._params[i]['TA'])
                                    qsptheta[i][j][kc] += exp(
                                        -0.5 * tmp * tmp)
                                    norms[i][j][kc] += 1
                            elif t in ["pent_plan", "pent_plan_max"]:
                                tmp = 0.4 if thetak <= self._params[i][
                                    'TA'] * pi \
                                    else 0.8
                                tmp2 = self._params[i]['IGW_TA'] * (
                                        thetak * ipi - tmp)
                                gaussthetak[i] = exp(-0.5 * tmp2 * tmp2)
                                if t == "pent_plan_max":
                                    qsptheta[i][j][kc] += gaussthetak[i]
                                    norms[i][j][kc] += 1
                            elif t == "bcc" and j < k:
                                if thetak >= self._params[i]['min_SPP']:
                                    tmp = self._params[i]['IGW_SPP'] * (
                                            thetak * ipi - 1.0)
                                    qsptheta[i][j][kc] += (
                                                self._params[i]['w_SPP'] *
                                                exp(-0.5 * tmp * tmp))
                                    norms[i][j][kc] += self._params[i]['w_SPP']
                            elif t == "sq_face_cap_trig_pris":
                                if thetak < self._params[i]['TA3']:
                                    tmp = self._params[i]['IGW_TA1'] * (
                                            thetak * ipi - self._params[i][
                                        'TA1'])
                                    qsptheta[i][j][kc] += exp(-0.5 * tmp * tmp)
                                    norms[i][j][kc] += 1

                        for m in range(nneigh):
                            if (m != j) and (m != k) and (not flag_xaxis):
                                tmp = max(
                                    -1.0, min(np.inner(zaxis, rijnorm[m]), 1.0))
                                thetam = acos(tmp)
                                xtwoaxistmp = gramschmidt(rijnorm[m], zaxis)
                                l = np.linalg.norm(xtwoaxistmp)
                                if l < very_small:
                                    flag_xtwoaxis = True
                                else:
                                    xtwoaxis = xtwoaxistmp / l
                                    phi = acos(max(
                                        -1.0,
                                        min(np.inner(xtwoaxis, xaxis), 1.0)))
                                    flag_xtwoaxis = False
                                    if self._comp_azi:
                                        phi2 = atan2(
                                            np.dot(xtwoaxis, yaxis),
                                            np.dot(xtwoaxis, xaxis))
                                        # print('{} {}'.format(180*phi/pi, 180*phi2/pi))
                                # South pole contributions of m.
                                if t in ["tri_bipyr", "sq_bipyr", "pent_bipyr",
                                         "hex_bipyr", "oct_max", "sq_plan_max",
                                         "hex_plan_max", "see_saw_rect"]:
                                    if thetam >= self._params[i]['min_SPP']:
                                        tmp = self._params[i]['IGW_SPP'] * (
                                                thetam * ipi - 1.0)
                                        qsptheta[i][j][kc] += exp(
                                            -0.5 * tmp * tmp)
                                        norms[i][j][kc] += 1

                                # Contributions of j-i-m angle and
                                # angles between plane j-i-k and i-m vector.
                                if not flag_xaxis and not flag_xtwoaxis:
                                    for i, t in enumerate(self._types):
                                        if t in ["tri_plan", "tri_plan_max",
                                                 "tet", "tet_max"]:
                                            tmp = self._params[i]['IGW_TA'] * (
                                                    thetam * ipi -
                                                    self._params[i]['TA'])
                                            tmp2 = cos(
                                                self._params[i]['fac_AA'] *
                                                phi) ** self._params[i][
                                                       'exp_cos_AA']
                                            tmp3 = 1 if t in ["tri_plan_max",
                                                              "tet_max"] \
                                                else gaussthetak[i]
                                            qsptheta[i][j][kc] += tmp3 * exp(
                                                -0.5 * tmp * tmp) * tmp2
                                            norms[i][j][kc] += 1
                                        elif t in ["pent_plan",
                                                   "pent_plan_max"]:
                                            tmp = 0.4 if thetam <= \
                                                         self._params[i][
                                                             'TA'] * pi \
                                                else 0.8
                                            tmp2 = self._params[i]['IGW_TA'] * (
                                                    thetam * ipi - tmp)
                                            tmp3 = cos(phi)
                                            tmp4 = 1 if t == "pent_plan_max" \
                                                else gaussthetak[i]
                                            qsptheta[i][j][kc] += tmp4 * exp(
                                                -0.5 * tmp2 * tmp2) * tmp3 * tmp3
                                            norms[i][j][kc] += 1
                                        elif t in ["T", "tri_pyr", "sq_pyr",
                                                   "pent_pyr", "hex_pyr"]:
                                            tmp = cos(
                                                self._params[i]['fac_AA'] *
                                                phi) ** self._params[i][
                                                      'exp_cos_AA']
                                            tmp3 = self._params[i]['IGW_EP'] * (
                                                    thetam * ipi - 0.5)
                                            qsptheta[i][j][kc] += tmp * exp(
                                                -0.5 * tmp3 * tmp3)
                                            norms[i][j][kc] += 1
                                        elif t in ["sq_plan", "oct",
                                                   "oct_legacy"]:
                                            if thetak < self._params[i][
                                                'min_SPP'] and \
                                                    thetam < self._params[i][
                                                'min_SPP']:
                                                tmp = cos(
                                                    self._params[i]['fac_AA'] *
                                                    phi) ** self._params[i][
                                                          'exp_cos_AA']
                                                tmp2 = self._params[i][
                                                           'IGW_EP'] * (
                                                               thetam * ipi - 0.5)
                                                qsptheta[i][j][kc] += tmp * exp(
                                                    -0.5 * tmp2 * tmp2)
                                                if t == "oct_legacy":
                                                    qsptheta[i][j][kc] -= tmp * \
                                                                          self._params[
                                                                              i][
                                                                              6] * \
                                                                          self._params[
                                                                              i][
                                                                              7]
                                                norms[i][j][kc] += 1
                                        elif t in ["tri_bipyr", "sq_bipyr",
                                                   "pent_bipyr",
                                                   "hex_bipyr", "oct_max",
                                                   "sq_plan_max",
                                                   "hex_plan_max"]:
                                            if thetam < self._params[i][
                                                'min_SPP']:
                                                if thetak < self._params[i][
                                                    'min_SPP']:
                                                    tmp = cos(self._params[i][
                                                                  'fac_AA'] *
                                                              phi) ** \
                                                          self._params[i][
                                                              'exp_cos_AA']
                                                    tmp2 = self._params[i][
                                                               'IGW_EP'] * (
                                                                   thetam * ipi - 0.5) if t != "hex_plan_max" else \
                                                        self._params[i][
                                                            'IGW_TA'] * (
                                                                fabs(
                                                                    thetam * ipi - 0.5) -
                                                                self._params[i][
                                                                    'TA'])
                                                    qsptheta[i][j][
                                                        kc] += tmp * exp(
                                                        -0.5 * tmp2 * tmp2)
                                                    norms[i][j][kc] += 1
                                        elif t == "bcc" and j < k:
                                            if thetak < self._params[i][
                                                'min_SPP']:
                                                if thetak > piover2:
                                                    fac = 1.0
                                                else:
                                                    fac = -1.0
                                                tmp = (thetam - piover2) / asin(
                                                    1 / 3)
                                                qsptheta[i][j][kc] += fac * cos(
                                                    3.0 * phi) * fac_bcc * \
                                                                      tmp * exp(
                                                    -0.5 * tmp * tmp)
                                                norms[i][j][kc] += 1
                                        elif t == "see_saw_rect":
                                            if thetam < self._params[i][
                                                'min_SPP']:
                                                if thetak < self._params[i][
                                                    'min_SPP'] and phi < 0.75 * pi:
                                                    tmp = cos(self._params[i][
                                                                  'fac_AA'] *
                                                              phi) ** \
                                                          self._params[i][
                                                              'exp_cos_AA']
                                                    tmp2 = self._params[i][
                                                               'IGW_EP'] * (
                                                                   thetam * ipi - 0.5)
                                                    qsptheta[i][j][kc] += tmp * \
                                                                          exp(
                                                                              -0.5 * tmp2 * tmp2)
                                                    norms[i][j][kc] += 1.0
                                        elif t in ["cuboct", "cuboct_max"]:
                                            if thetam < self._params[i][
                                                'min_SPP'] and \
                                                    thetak > self._params[i][
                                                4] and \
                                                    thetak < self._params[i][2]:
                                                if thetam > self._params[i][
                                                    4] and \
                                                        thetam < \
                                                        self._params[i][2]:
                                                    tmp = cos(phi)
                                                    tmp2 = self._params[i][
                                                               5] * (
                                                                       thetam * ipi - 0.5)
                                                    qsptheta[i][j][
                                                        kc] += tmp * tmp * exp(
                                                        -0.5 * tmp2 * tmp2)
                                                    norms[i][j][kc] += 1.0
                                                elif thetam < self._params[i][
                                                    4]:
                                                    tmp = 0.0556 * (cos(
                                                        phi - 0.5 * pi) - 0.81649658)
                                                    tmp2 = self._params[i][
                                                               6] * (
                                                                   thetam * ipi - onethird)
                                                    qsptheta[i][j][kc] += exp(
                                                        -0.5 * tmp * tmp) * \
                                                                          exp(
                                                                              -0.5 * tmp2 * tmp2)
                                                    norms[i][j][kc] += 1.0
                                                elif thetam > self._params[i][
                                                    2]:
                                                    tmp = 0.0556 * (cos(
                                                        phi - 0.5 * pi) - \
                                                                    0.81649658)
                                                    tmp2 = self._params[i][
                                                               6] * (
                                                                       thetam * ipi - \
                                                                       twothird)
                                                    qsptheta[i][j][kc] += exp(
                                                        -0.5 * tmp * tmp) * \
                                                                          exp(
                                                                              -0.5 * tmp2 * tmp2)
                                                    norms[i][j][kc] += 1.0
                                        elif t == "sq_face_cap_trig_pris" and not flag_yaxis:
                                            if thetak < self._params[i]['TA3']:
                                                if thetam < self._params[i][
                                                    'TA3']:
                                                    tmp = cos(self._params[i][
                                                                  'fac_AA1'] * \
                                                              phi2) ** \
                                                          self._params[i][
                                                              'exp_cos_AA1']
                                                    tmp2 = self._params[i][
                                                               'IGW_TA1'] * (
                                                                   thetam * ipi -
                                                                   self._params[
                                                                       i][
                                                                       'TA1'])
                                                else:
                                                    tmp = cos(self._params[i][
                                                                  'fac_AA2'] * \
                                                              (phi2 +
                                                               self._params[i][
                                                                   'shift_AA2'])) ** \
                                                          self._params[i][
                                                              'exp_cos_AA2']
                                                    tmp2 = self._params[i][
                                                               'IGW_TA2'] * (
                                                                   thetam * ipi -
                                                                   self._params[
                                                                       i][
                                                                       'TA2'])
                                                    # print("phi2 {}   phi2+shift {}   tmp {}  ".format(phi2*180/pi, 180*(phi2 + self._params[i]['shift_AA2'])/pi, tmp))
                                                qsptheta[i][j][kc] += tmp * exp(
                                                    -0.5 * tmp2 * tmp2)
                                                norms[i][j][kc] += 1

                        kc += 1

            # Normalize Peters-style OPs.
            for i, t in enumerate(self._types):
                if t in ["tri_plan", "tet", "bent", "sq_plan",
                         "oct", "oct_legacy", "cuboct", "pent_plan"]:
                    ops[i] = tmp_norm = 0.0
                    for j in range(nneigh):
                        ops[i] += sum(qsptheta[i][j])
                        tmp_norm += float(sum(norms[i][j]))
                    ops[i] = ops[i] / tmp_norm if tmp_norm > 1.0e-12 else None
                elif t in ["T", "tri_pyr", "see_saw_rect", "sq_pyr",
                           "tri_bipyr",
                           "sq_bipyr", "pent_pyr", "hex_pyr", "pent_bipyr",
                           "hex_bipyr", "oct_max", "tri_plan_max", "tet_max",
                           "sq_plan_max", "pent_plan_max", "cuboct_max",
                           "hex_plan_max",
                           "sq_face_cap_trig_pris"]:
                    ops[i] = None
                    if nneigh > 1:
                        for j in range(nneigh):
                            for k in range(len(qsptheta[i][j])):
                                qsptheta[i][j][k] = qsptheta[i][j][k] / \
                                                    norms[i][j][k] \
                                    if norms[i][j][k] > 1.0e-12 else 0.0
                            ops[i] = max(qsptheta[i][j]) if j == 0 \
                                else max(ops[i], max(qsptheta[i][j]))
                elif t == "bcc":
                    ops[i] = 0.0
                    for j in range(nneigh):
                        ops[i] += sum(qsptheta[i][j])
                    ops[i] = ops[i] / float(0.5 * float(
                        nneigh * (6 + (nneigh - 2) * (nneigh - 3)))) \
                        if nneigh > 3 else None
                elif t == "sq_pyr_legacy":
                    if nneigh > 1:
                        dmean = np.mean(dist)
                        acc = 0.0
                        for d in dist:
                            tmp = self._params[i][2] * (d - dmean)
                            acc = acc + exp(-0.5 * tmp * tmp)
                        for j in range(nneigh):
                            ops[i] = max(qsptheta[i][j]) if j == 0 \
                                else max(ops[i], max(qsptheta[i][j]))
                        ops[i] = acc * ops[i] / float(nneigh)
                        # nneigh * (nneigh - 1))
                    else:
                        ops[i] = None

        # Then, deal with the new-style OPs that require vectors between
        # neighbors.
        if self._geomops2:
            # Compute all (unique) angles and sort the resulting list.
            aij = []
            for ir, r in enumerate(rijnorm):
                for j in range(ir + 1, len(rijnorm)):
                    aij.append(
                        acos(max(-1.0, min(np.inner(r, rijnorm[j]), 1.0))))
            aijs = sorted(aij)

            # Compute height, side and diagonal length estimates.
            neighscent = np.array([0.0, 0.0, 0.0])
            for j, neigh in enumerate(neighsites):
                neighscent = neighscent + neigh.coords
            if nneigh > 0:
                neighscent = (neighscent / float(nneigh))
            h = np.linalg.norm(neighscent - centvec)
            b = min(distjk_unique) if len(distjk_unique) > 0 else 0
            dhalf = max(distjk_unique) / 2.0 if len(distjk_unique) > 0 else 0

            for i, t in enumerate(self._types):
                if t == "reg_tri" or t == "sq":
                    if nneigh < 3:
                        ops[i] = None
                    else:
                        ops[i] = 1.0
                        if t == "reg_tri":
                            a = 2.0 * asin(b / (2.0 * sqrt(h * h + (b / (
                                    2.0 * cos(3.0 * pi / 18.0))) ** 2.0)))
                            nmax = 3
                        elif t == "sq":
                            a = 2.0 * asin(
                                b / (2.0 * sqrt(h * h + dhalf * dhalf)))
                            nmax = 4
                        for j in range(min([nneigh, nmax])):
                            ops[i] = ops[i] * exp(-0.5 * ((
                                                                  aijs[j] - a) *
                                                          self._params[i][
                                                              0]) ** 2)

        return ops