Periodic Bloch-Floquet boundary conditions for calculating phononic band stucture

I am trying to calculate the phononic band structure of a thin cross-hole perforated silicon phononic crystal using the partial differential equation toolbox. However, I don't know how to set-up periodic Bloch-Floquet boundary conditions so that I can loop over a set of kx and ky wavevectors. Sofar my implementation looks as shown below, where I calculate eigenmodes for the unit cell with outer edges held constant at u=0. Thank you.

a=500e-9 %box size/ lattice pitch
E=170e9 %Youngs modulus in [Pa]
nu=0.28 %Poission's ratio
rho=2329 %mass density in [kg/m^3]

x_sqr=[ 0.5*a -0.5*a -0.5*a 0.5*a]; y_sqr=[-0.5*a -0.5*a  0.5*a 0.5*a];
sqr=[3 4 x_sqr y_sqr]';

x_R1=[0.4*a -0.4*a -0.4*a 0.4*a]; y_R1=[-0.1*a -0.1*a  0.1*a 0.1*a];
R1=[3 4 x_R1 y_R1]';

x_R2=[-0.1*a -0.1*a  0.1*a 0.1*a]; y_R2=[0.4*a -0.4*a -0.4*a 0.4*a];
R2=[3 4 x_R2 y_R2]'; 

ns = char('sqr','R1','R2'); ns=ns';
gd = [sqr,R1,R2];
sf = '(sqr-(R1+R2))';
[geometry,bool]=decsg(gd,sf,ns);
outer_edges=[1 2 14 15]'

N = 2;
model = createpde(N);

geometryFromEdges(model,geometry)
generateMesh(model,'Hmax',0.04*a);
applyBoundaryCondition(model,'Edge',outer_edges,'u',0,'EquationIndex',2)

figure;
pdegplot(model,'EdgeLabels','on');

figure;
pdeplot(model)

%Plane strain c matrix
G=E/(2*(1+nu))
mu=2*G*(nu/(1-nu))
c11=[2*G+mu 0 G];
c12=[0 mu G 0];
c22=[G 0 2*G+mu];
C=[c11 c12 c22]';

A = 0
D = rho;
R = [-0.1 2e+21]

[eVec,eVal] = pdeeig(model,C,A,D,R); % Solve PDE
numEig = size(eVal,1);
u2 = reshape(eVec,[],2,numEig);

for i=1:numEig
 figure('Name',['Eigenmode ' num2str(i)]);
 pdeplot(model,'xydata',u2(:,2,i),'zdata',u2(:,2,i))
 axis tight
end