Defining Node Sets of Physical Groups using Python API

I am creating a simple finite element model of a cylinder with spherical holes being cut out uniformly throughout the model. This model is meant to be analyzed using ABAQUS. So, I defined the volume itself and the top, bottom, and middle surfaces of the cylinder using physical groups and then used "gmsh.option.setNumber('Mesh.SaveGroupsOfNodes',1)" to save nodes in each physical group in the outputted INP file. The problem I'm having is that the nodes that are being outputted onto the INP file are not correct. When I run this model without the holes, the correct nodes are outputted, but incorrect nodes are given when the holes are added in.

I believe there is a simple fix that I am missing, but I am unable to find this fix. Enclosed below is a link to the code I made using the Python API.

holes = []
x = 0
y = 0
z = 0
tot = 0
theta = 0
z_pos = 0.75

import gmsh
import math
import sys

gmsh.initialize()
gmsh.logger.start()
gmsh.model.add("Figure 2D GMSH Model - Uniform Cylinder")
gmsh.model.occ.addCylinder(0,0,0,0,0,15,6.5,1)
gmsh.model.occ.synchronize()

while z_pos <= 15:
    # Center hole
    tot += 1
    gmsh.model.occ.addSphere(0,0,z_pos,0.75,1+tot)
    holes.append((3,1+tot))
    theta = 0

    for hole in range(6):
        tot += 1
        x = 1.625 * math.cos(math.radians(theta))
        y = 1.625 * math.sin(math.radians(theta))
        gmsh.model.occ.addSphere(x, y, z_pos, 0.75, 1 + tot)
        holes.append((3, 1 + tot))
        theta += 60

    theta = 0
    for hole in range(12):
        tot += 1
        x = 3.25 * math.cos(math.radians(theta))
        y = 3.25 * math.sin(math.radians(theta))
        gmsh.model.occ.addSphere(x, y, z_pos, 0.75, 1 + tot)
        holes.append((3, 1 + tot))
        theta += 30

    theta = 0
    for hole in range(20):
        tot += 1
        x = 4.875 * math.cos(math.radians(theta))
        y = 4.875 * math.sin(math.radians(theta))
        gmsh.model.occ.addSphere(x, y, z_pos, 0.75, 1 + tot)
        holes.append((3, 1 + tot))
        theta += 18

    theta = 0
    for hole in range(24):
        tot += 1
        x = 6.5 * math.cos(math.radians(theta))
        y = 6.5 * math.sin(math.radians(theta))
        gmsh.model.occ.addSphere(x, y, z_pos, 0.75, 1 + tot)
        holes.append((3, 1 + tot))
        theta += 15
    gmsh.model.occ.synchronize()
    z_pos += 1.5

gmsh.model.addPhysicalGroup(3, [1], 1000)
gmsh.model.addPhysicalGroup(2, [1], 1001)
gmsh.model.addPhysicalGroup(2, [2], 1002)
gmsh.model.addPhysicalGroup(2, [3], 1003)
gmsh.model.setPhysicalName(3,1000,"cube")
gmsh.model.setPhysicalName(2,1001,"middle_surface")
gmsh.model.setPhysicalName(2,1002,"top_surface")
gmsh.model.setPhysicalName(2,1003,"bottom_surface")
gmsh.option.setNumber('Mesh.SaveGroupsOfNodes',1)

gmsh.model.mesh.setSize(gmsh.model.getEntities(0), 1)

gmsh.model.mesh.generate(3)
gmsh.write("Figure 2D GMSH Model Duplicate- newattemp2.inp")
# Inspect the log:
log = gmsh.logger.get()
print("Logger has recorded " + str(len(log)) + " lines")
gmsh.logger.stop()
if '-nopopup' not in sys.argv:
    gmsh.fltk.run()
gmsh.finalize()

Kindly, Bobsy Narayan