def sqdist(point1, point2, rotmat):
    """
    This routine calculates the anisotropic distance between two points
    given the coordinates of each point and a definition of the
    anisotropy.

    This method only consider a single anisotropy senario.

    Parameters
    ----------
    point1 : tuple
        Coordinates of first point (x1,y1,z1)
    point2 : tuple
        Coordinates of second point (x2,y2,z2)
    rotmat : 3*3 ndarray
        matrix of rotation for this structure

    Returns
    -------
    sqdist : scalar
        The squared distance accounting for the anisotropy
        and the rotation of coordinates (if any).
    """
    dx = point1[0] - point2[0]
    dy = point1[1] - point2[1]
    dz = point1[2] - point2[2]
    sqdist = 0.0
    for i in range(3):
        cont = rotmat[i, 0] * dx + \
                rotmat[i, 1] * dy + \
                rotmat[i, 2] * dz
        sqdist += cont * cont
    return sqdist