The long-range corrected spin-unrestricted density functional theory, LC-UBLYP, method has been employed to unravel the relationship between the hole (referred to as "antidot") structure of hexagonal graphene nanoflakes (HGNFs) and their properties, including their multiradical characters y [the occupation number of the lowest unoccupied natural orbital (LUNO)+i (i = 0, 1,...)], aromaticity, and the second hyperpolarizabilities (γ), which is the third-order nonlinear optical (NLO) response at the molecular level. These systems exhibit a wide range of open-shell/multiradical character, which shows an oscillatory behavior as a function of increasing the size of the antidot. This structural dependence is shown to originate from the oscillatory variations in the HOMO-LUMO energy gap, which result from the fact that the bonding and antibonding interactions between the central antidot-shape HGNFs and the surroundings alternate in both the HOMO and the LUMO as going from the center of the molecule to the peripheral region. Moreover, the multiradical character of these structures is strongly correlated with the variations of magnetic criteria of aromaticity, the nucleus-independent chemical shift (NICS) and the out-of-plane diagonal component of the magnetic shielding tensor (-σ ), as well as with the γ. So, systems with larger multiradical characters tend to exhibit larger NICS and - σ values, i.e., less aromaticity. Although the number of π electrons of the antidot HGNFs is smaller than that of the corresponding perfect C H HGNF, their γ value can be larger as evidenced by the antidot C H HGNF, which has intermediate multiradical characters (y0 = y1 = 0.658) and displays a γvalue 1.7 times larger than that of the perfect HGNF. These strong impacts of the antidot structure on the third-order NLO properties indicate that the antidot HGNFs are promising building blocks for a new class of multiradical NLO materials, the properties of which can be controlled by adjusting the antidot size.