creating hexagonal lattice graph in networkx

I need to create a hexagonal lattice network based on another researcher's work. The researcher describes describes the structure of their network in the yellow highlighted section of text here.

Here's a little bit more info from the actual paper (again, I'm only in need of help with the hexagonal lattice construction):

The full paper is here.

I suspect I'm missing something. Here's my code:

import networkx as nx
import random
import matplotlib.pyplot as plt

def create_hexagonal_lattice_graph(rows, cols):
    G = nx.hexagonal_lattice_graph(rows, cols)
    simple_G = nx.Graph(G)
    return simple_G

# Define the number of rows and columns to get approximately 200 nodes
rows, cols = 6, 6  # Adjusting to get roughly 200 nodes (10*10 grid)

# Create the initial hexagonal lattice graph
G = create_hexagonal_lattice_graph(rows, cols)

# Ensure all nodes have degree 6 by adding missing edges
for node in G.nodes():
    while G.degree(node) < 6:
        possible_nodes = set(G.nodes()) - set(G.neighbors(node)) - {node}
        new_neighbor = random.choice(list(possible_nodes))
        G.add_edge(node, new_neighbor)

# Check the number of nodes and ensure all have degree 6
num_nodes = G.number_of_nodes()
degrees = dict(G.degree())
degree_check = all(degree == 6 for degree in degrees.values())

print(f"Number of nodes: {num_nodes}")
print(f"All nodes have degree 6: {degree_check}")

plt.title("Initial Hexagonal Lattice")
nx.draw_circular(G, node_size=20, with_labels=False)

The result does not have the degree distribution and clustering coefficient that I am hoping for. I feel like I am missing necessary information to construct the network. How do I modify my code to create an undirected hexagonal lattice graph that satisfies these bare-bones characteristics (Z=6; CC=.4)?

IIUC, you could approach it geometrically by using a DistanceBand with Chebyshev's distance :

import networkx as nx
from libpysal import weights

def lattice(r, c, no_ortho):
    """Return Graph and nodes' positions."""

    from itertools import product
    grid = list(product(range(r), range(c)))
    if no_ortho: grid = grid[1::2]
        
    G = weights.DistanceBand(
        grid,
        threshold=int(no_ortho) + 1,
        p=float("inf"),  # chebyshev
        silence_warnings=True,
    ).to_networkx()
    
    if no_ortho:
        G.remove_edges_from(
            [(u, v) for u, v in G.edges if abs(u - v) == 1]
        )
        
    return G, dict(zip(G, grid))

A = lattice(7, 7, True)
B = lattice(7, 7, False)