The paper you referenced has detailed documentation on the authors' simulation methodologies- I would suggest you re-read that more carefully as there are several discrepancies between your simulation and theirs. If results still differ once the methodologies are identical, then there is something to look into.
Most notably, if you're trying to observe out-of-plane buckling/wrinkles, how will the graphene sheet accommodate them if you've fixed the z dimension?
I understand that my approach is different from the one described in the paper which I have referenced but nonetheless my results should still reflect those obtained by them, that is a decreasing sheet size with increasing percentage defects.
So, to put it another way, how may I be sure that a graphene sheet with certain atoms missing is properly minimized if I am getting practically the same sheet size irrespective of how many atoms I’m deleting?
A note of clarification: What I actually meant was that the z-direction has a fixed boundary condition applied to it in contrast to the x and y dimensions which are periodic. This was done in an attempt to simulate bulk graphene. The length of the z-dimension is 100 Å which should be more than enough for the system to buckle.
Sorry for the misunderstanding.