Exact diagonalization of graphene's tight binding Hamiltonian While directly diagonalize graphene's tight binding Hamiltonian, which is numerical. We have to use a finite-sized graphene.
So how to deal with boundary conditions? The usual solutions are zigzag or armchair condition, but to make our model more realistic to real infinite graphene plane, how about using periodic condition at boundaries while exactly diagonalizing the tight binding Hamiltonian? 
 A: You're mixing two things here. One is what the structure of the boundary is, e.g. armchair or zigzag. The other is what the wavefunction does at the boundary.
For your finite size cluster of carbon atoms, you have to decide what shape it has, which basically means deciding how many lattice sites you include and where you put them. This would decide whether your finite cluster has armchair or zigzag boundaries.
If  you treat these boundaries as "real" system boundaries, then you have closed boundary condition. But if you want periodic boundary conditions, you have to decide on a way to map one boundary to the other, i.e., how does everything get "wrapped around". 
E.g. for a rectangular lattice this is very easy, you just declare that electrons leaving the lattice to the left re-enter it on the right, basically stating that site $0$ and site $N$ are the same.
For your graphene lattice, you can do the same, you just have to be careful that the boundary of your finite sized cluster is such that this wrapping actually works. That means: For every atom on the boundary of your finite cell, there must still be a way to identify its 6 unique "nearest neighbors", some of which are "actual" neighbors, and some appear on the opposite site of the cluster.
A: TL;DR version: you want periodic boundary conditions with an extra twist.
You don't want straightforward periodic boundary conditions as this would be solving for a structure that is periodic. You can solve for an infinite sheet of graphene using Bloch's theorem.
I'll give a couple of details:
Find a basic finite unit that you can tile to make the infinite sheet so that you cover the plane by placing copies of the tile at positions $ma_1+na_2$ with one copy for each integer pair $(m,n)$.
Using Bloch's theorem, a basis for the space of wavefunctions can be found as follows:
Pick phase angles $k_1$ and $k_2$. Now imagine we're working with periodic boundary conditions so the hopping matrix has elements corresponding to neighbouring pairs of atoms where the elements of the pair are on opposite sides of the tile. When you walk off one side and come back on the other, this must correspond to a step along a vector of the form $ma_1+na_2$ for some $m$ and $n$. You modify the hopping matrix so that it has an extra phase factor $\exp(mk_1+nk_2)$ (and the complex conjugate for the step going the other way). You now solve the tight-binding model numerically using this new "twisted" matrix. As you're dealing with a finite tile, if it has $n$ atoms then you expect $n$ eigenfunctions.
You get a basis for the whole plane by considering all possible wavefunctions as you vary $k_1$ and $k_2$, ie. you get $n$ bands with each band parameterised by $k_1$ and $k_2$.
For graphene you can use a repeating two atom unit: one atom has a bond going straight out to the left and the other has one going straight out to the right.
Have a look at section 2.1.2 in this course.
