Generate a random graph with geometrical degree distribution
You can use the configuration model for this. Create a "stub" array where each of the $n$ vertices appears in the array $k_i \sim \text{Geom}(p)$ times. Shuffle the array and then treat each adjacent pair in the array as an edge in the graph. Here's some python code to illustrate.

from numpy.random.mtrand import geometric
from random import shuffle
from _collections import defaultdict

if __name__ == '__main__':
    # Number of vertices
    n = 1000

    # Geometric distribution probability parameter
    p = 0.5

    stubs = []
    for i in xrange(0, n):
        # Sample the degree from the geometric distribution
        degree = geometric(p) 
        # Insert into the stubs array the current vertex label degree times
        stubs += [i] * degree

    # Shuffle the values
    shuffle(stubs)

    # Build the graph
    G = {}
    for i in xrange(0, n):
        G[i] = []

    for i in xrange(1, len(stubs), 2):
        G[stubs[i]].append(stubs[i-1])
        G[stubs[i-1]].append(stubs[i])

    # Recover the degree distribution
    degrees = { i : len(neighbors) for i, neighbors in G.items() }

    dist = defaultdict(int)
    for (vertex, degree) in degrees.items():
        dist[degree] += 1

    print("Degree, empirical, theoretical")
    for (degree, freq) in dist.items():
        print("%d %f %f" % (degree, freq / float(n),  pow((1 - p), degree - 1) * p ))


Sample output:

Degree, empirical, theoretical
1 0.496000 0.500000
2 0.254000 0.250000
3 0.128000 0.125000
4 0.071000 0.062500
5 0.023000 0.031250
6 0.018000 0.015625
7 0.005000 0.007812
8 0.003000 0.003906
9 0.002000 0.001953


For more detailed description of the the configuration model, see Section 13.2 p. 435 of Networks by Newman.