def generateMonkhorstPackGrid (n_abc, R, time_reversal_symmetry):
    import math
    from numpy import array, arange, ones

    '''
    Returns the k-points of a (na x nb x nc) Monkhorst-Pack grid
    for (reciprocal) lattice vectors R = [A,B,C]
    where na = n_abc[0] etc
    Optionally, the grid can be reduced by time-reversal symmetry,
    which removes -k where k is present, except for the Gamma point
    The general formulae for the grid are Eqs. 3-4 in
    H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976)
    which state that
        (3)  ur = (2r-N-1)/2N   where   r=1,2,...,N
    where N is the number of states in each direction.
    These integers are passed as na, nb, nc to this function
    The k-points are then
        (4)  k = up * A + uq * B + ur * C
    where A, B, C are the reciprocal lattice vectors
    '''
    
    '''
    In the case of time-reversal symmetry, there are half as many
    k-points plus the Gamma point, if included. The total number
    of points na*nb*nc is odd if the Gamma point is included,
    thus (na*nb*nc+1)/2 using integer division will give the
    correct number of points.
    '''
    n_product = n_abc[0]*n_abc[1]*n_abc[2]
    total_number_of_kpoints = n_product
    if (time_reversal_symmetry):
        total_number_of_kpoints = (n_product+1)/2

    # Create the numbers ur according to Eq. (3)
    ua = (2*arange(1,n_abc[0]+1)-n_abc[0]-1)/(2.*n_abc[0])
    ub = (2*arange(1,n_abc[1]+1)-n_abc[1]-1)/(2.*n_abc[1])
    uc = (2*arange(1,n_abc[2]+1)-n_abc[2]-1)/(2.*n_abc[2])

    if (time_reversal_symmetry):
        '''
        Impose time-reversal symmetry by removing one half of ua
        (or half of ub, if ua has only a single point, 
        or half of uc if both ua and ub have a single point, 
        or remove none if there is just a single point)
        '''
        if   len(ua)>1:
            ua = ua[0:(len(ua)+1)/2]
        elif len(ub)>1:
            ub = ub[0:(len(ub)+1)/2]
        elif len(uc)>1:
            uc = uc[0:(len(uc)+1)/2]

    # Weights
    if (time_reversal_symmetry):
        wG = 1./n_product
        w0 = 2.*wG
        # If na*nb*nc is odd, it means that the Gamma
        # point (0,0,0) is included.
        # It has half of the weight of the other points
        # if time-reversal symmetry is used
        if (n_product % 2):
            W = [w0,]*(total_number_of_kpoints-1)+[wG]
        else:
            W = [w0,]*total_number_of_kpoints
    else:
        w0 = 1./(total_number_of_kpoints)
        W = [w0,]*total_number_of_kpoints
    
    # K-points
    K = []
    for i in range(len(ua)):
        for j in range(len(ub)):
            for k in range(len(uc)):
                # Calculate k-points according to Eq. (4)
                K.append(ua[i]*R[0]+ub[j]*R[1]+uc[k]*R[2])
                '''
                When time-reversal is used, the ua vectors are ordered
                such that when we hit the Gamma point, the remaining
                points are time-reversed copies of the already
                included points, and should not be included
                '''
                if (time_reversal_symmetry) & \
                       (ua[i]==0) & (ub[j]==0) & (uc[k]==0):
                    return K,W
    return K,W

def printMonkhorstPackGrid(K,W):
    for i in range(len(K)):
        print "%.5f  \t%.5f  \t%.5f  \t%.5f" %\
              (K[i][0],K[i][1],K[i][2],W[i])

def MonkhorstPackGrid (n_abc,
                       bulk_configuration,
                       time_reversal_symmetry=True,
                       verbosity_level=0):
    from bulk_utilities import reciprocalVectors
    R = reciprocalVectors(bulk_configuration)
    (K,W) = generateMonkhorstPackGrid (n_abc, R, time_reversal_symmetry)
    if verbosity_level>=10:
        printMonkhorstPackGrid(K,W)
    return (K,W)