diff --git a/openmc/surface.py b/openmc/surface.py
index 38a2c92d84..6ce0fe0da5 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -1,6 +1,8 @@
+from __future__ import division
 from abc import ABCMeta
 from collections import Iterable, OrderedDict
 from copy import deepcopy
+from functools import partial
 from numbers import Real, Integral
 from xml.etree import ElementTree as ET
 from math import sqrt
@@ -1326,6 +1328,7 @@ class Sphere(Surface):
         z = point[2] - self.z0
         return x**2 + y**2 + z**2 - self.r**2
 
+
 @add_metaclass(ABCMeta)
 class Cone(Surface):
     """A conical surface parallel to the x-, y-, or z-axis.
@@ -1995,147 +1998,66 @@ def get_rectangular_prism(width, height, axis='z', origin=(0., 0.),
     check_type('width', width, Real)
     check_type('height', height, Real)
     check_type('corner_radius', corner_radius, Real)
-    check_value('axis', axis, ['x','y','z'])
+    check_value('axis', axis, ['x', 'y', 'z'])
     check_type('origin', origin, Iterable, Real)
 
+    # Define function to create a plane on given axis
+    def plane(axis, name, value):
+        cls = globals()['{}Plane'.format(axis.upper())]
+        return cls(name='{} {}'.format(name, axis),
+                   boundary_type=boundary_type,
+                   **{axis + '0': value})
+
     if axis == 'x':
-        min_y = YPlane(name='minimum y', y0=-width/2.+origin[0],
-                       boundary_type=boundary_type)
-        max_y = YPlane(name='maximum y', y0=+width/2.+origin[0],
-                       boundary_type=boundary_type)
-        min_z = ZPlane(name='minimum z', z0=-height/2.+origin[1],
-                       boundary_type=boundary_type)
-        max_z = ZPlane(name='maximum z', z0=+height/2.+origin[1],
-                       boundary_type=boundary_type)
-        prism = +min_y & -max_y & +min_z & -max_z
+        x1, x2 = 'y', 'z'
     elif axis == 'y':
-        min_x = XPlane(name='minimum x', x0=-width/2.+origin[0],
-                       boundary_type=boundary_type)
-        max_x = XPlane(name='maximum x', x0=+width/2.+origin[0],
-                       boundary_type=boundary_type)
-        min_z = ZPlane(name='minimum z', z0=-height/2.+origin[1],
-                       boundary_type=boundary_type)
-        max_z = ZPlane(name='maximum z', z0=+height/2.+origin[1],
-                       boundary_type=boundary_type)
-        prism = +min_x & -max_x & +min_z & -max_z
+        x1, x2 = 'x', 'z'
     else:
-        min_x = XPlane(name='minimum x', x0=-width/2.+origin[0],
-                       boundary_type=boundary_type)
-        max_x = XPlane(name='maximum x', x0=+width/2.+origin[0],
-                       boundary_type=boundary_type)
-        min_y = YPlane(name='minimum y', y0=-height/2.+origin[1],
-                       boundary_type=boundary_type)
-        max_y = YPlane(name='maximum y', y0=+height/2.+origin[1],
-                       boundary_type=boundary_type)
-        prism = +min_x & -max_x & +min_y & -max_y
+        x1, x2 = 'x', 'y'
 
+    # Get cylinder class corresponding to given axis
+    cyl = globals()['{}Cylinder'.format(axis.upper())]
+
+    # Create rectangular region
+    min_x1 = plane(x1, 'minimum', -width/2 + origin[0])
+    max_x1 = plane(x1, 'maximum', width/2 + origin[0])
+    min_x2 = plane(x2, 'minimum', -height/2 + origin[1])
+    max_x2 = plane(x2, 'maximum', height/2 + origin[1])
+    prism = +min_x1 & -max_x1 & +min_x2 & -max_x2
+
+    # Handle rounded corners if given
     if corner_radius > 0.:
-        if axis == 'x':
-            y_min_z_min = XCylinder(name='y min z min',
-                                    y0=origin[0] - width /2. + corner_radius,
-                                    z0=origin[1] - height/2. + corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
-            y_min_z_max = XCylinder(name='y min z max',
-                                    y0=origin[0] - width /2. + corner_radius,
-                                    z0=origin[1] + height/2. - corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
-            y_max_z_min = XCylinder(name='y max z min',
-                                    y0=origin[0] + width /2. - corner_radius,
-                                    z0=origin[1] - height/2. + corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
-            y_max_z_max = XCylinder(name='y max z max',
-                                    y0=origin[0] + width /2. - corner_radius,
-                                    z0=origin[1] + height/2. - corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
-            y_min = YPlane(name='min y', y0=-width/2.+origin[0]+corner_radius,
-                           boundary_type=boundary_type)
-            y_max = YPlane(name='max y', y0=+width/2.+origin[0]-corner_radius,
-                           boundary_type=boundary_type)
-            z_min = ZPlane(name='min z', z0=-height/2.+origin[1]+corner_radius,
-                           boundary_type=boundary_type)
-            z_max = ZPlane(name='max z', z0=+height/2.+origin[1]-corner_radius,
-                           boundary_type=boundary_type)
+        args = {'R': corner_radius, 'boundary_type': boundary_type}
 
-            corners = (+y_min_z_min & -y_min & -z_min) | \
-                      (+y_min_z_max & -y_min & +z_max) | \
-                      (+y_max_z_min & +y_max & -z_min) | \
-                      (+y_max_z_max & +y_max & +z_max)
+        args[x1 + '0'] = origin[0] - width/2 + corner_radius
+        args[x2 + '0'] = origin[1] - height/2 + corner_radius
+        x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args)
 
-        elif axis == 'y':
-            x_min_z_min = YCylinder(name='x min z min',
-                                    x0=origin[0] - width /2. + corner_radius,
-                                    z0=origin[1] - height/2. + corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
-            x_min_z_max = YCylinder(name='x min z max',
-                                    x0=origin[0] - width /2. + corner_radius,
-                                    z0=origin[1] + height/2. - corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
-            x_max_z_min = YCylinder(name='x max z min',
-                                    x0=origin[0] + width /2. - corner_radius,
-                                    z0=origin[1] - height/2. + corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
-            x_max_z_max = YCylinder(name='x max z max',
-                                    x0=origin[0] + width /2. - corner_radius,
-                                    z0=origin[1] + height/2. - corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
+        args[x1 + '0'] = origin[0] - width/2 + corner_radius
+        args[x2 + '0'] = origin[1] - height/2 + corner_radius
+        x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args)
 
-            x_min = XPlane(name='min x', x0=-width/2.+origin[0]+corner_radius,
-                           boundary_type=boundary_type)
-            x_max = XPlane(name='max x', x0=+width/2.+origin[0]-corner_radius,
-                           boundary_type=boundary_type)
-            z_min = ZPlane(name='min z', z0=-height/2.+origin[1]+corner_radius,
-                           boundary_type=boundary_type)
-            z_max = ZPlane(name='max z', z0=+height/2.+origin[1]-corner_radius,
-                           boundary_type=boundary_type)
+        args[x1 + '0'] = origin[0] - width/2 + corner_radius
+        args[x2 + '0'] = origin[1] + height/2 - corner_radius
+        x1_min_x2_max = cyl(name='{} min {} max'.format(x1, x2), **args)
 
-            corners = (+x_min_z_min & -x_min & -z_min) | \
-                      (+x_min_z_max & -x_min & +z_max) | \
-                      (+x_max_z_min & +x_max & -z_min) | \
-                      (+x_max_z_max & +x_max & +z_max)
+        args[x1 + '0'] = origin[0] + width/2 - corner_radius
+        args[x2 + '0'] = origin[1] - height/2 + corner_radius
+        x1_max_x2_min = cyl(name='{} max {} min'.format(x1, x2), **args)
 
-        else:
-            x_min_y_min = ZCylinder(name='x min y min',
-                                    x0=origin[0] - width /2. + corner_radius,
-                                    y0=origin[1] - height/2. + corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
-            x_min_y_max = ZCylinder(name='x min y max',
-                                    x0=origin[0] - width /2. + corner_radius,
-                                    y0=origin[1] + height/2. - corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
-            x_max_y_min = ZCylinder(name='x max y min',
-                                    x0=origin[0] + width /2. - corner_radius,
-                                    y0=origin[1] - height/2. + corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
-            x_max_y_max = ZCylinder(name='x max y max',
-                                    x0=origin[0] + width /2. - corner_radius,
-                                    y0=origin[1] + height/2. - corner_radius,
-                                    R=corner_radius,
-                                    boundary_type=boundary_type)
+        args[x1 + '0'] = origin[0] + width/2 - corner_radius
+        args[x2 + '0'] = origin[1] + height/2 - corner_radius
+        x1_max_x2_max = cyl(name='{} max {} max'.format(x1, x2), **args)
 
-            x_min = XPlane(name='min x', x0=-width/2.+origin[0]+corner_radius,
-                           boundary_type=boundary_type)
-            x_max = XPlane(name='max x', x0=+width/2.+origin[0]-corner_radius,
-                           boundary_type=boundary_type)
-            y_min = YPlane(name='min y', y0=-height/2.+origin[1]+corner_radius,
-                           boundary_type=boundary_type)
-            y_max = YPlane(name='max y', y0=+height/2.+origin[1]-corner_radius,
-                           boundary_type=boundary_type)
+        x1_min = plane(x1, 'min', -width/2 + origin[0] + corner_radius)
+        x1_max = plane(x1, 'max', width/2 + origin[0] - corner_radius)
+        x2_min = plane(x2, 'min', -height/2 + origin[1] + corner_radius)
+        x2_max = plane(x2, 'max', height/2 + origin[1] - corner_radius)
 
-            corners = (+x_min_y_min & -x_min & -y_min) | \
-                      (+x_min_y_max & -x_min & +y_max) | \
-                      (+x_max_y_min & +x_max & -y_min) | \
-                      (+x_max_y_max & +x_max & +y_max)
+        corners = (+x1_min_x2_min & -x1_min & -x2_min) | \
+                  (+x1_min_x2_max & -x1_min & +x2_max) | \
+                  (+x1_max_x2_min & +x1_max & -x2_min) | \
+                  (+x1_max_x2_max & +x1_max & +x2_max)
 
         prism = prism & ~corners
 
@@ -2170,8 +2092,8 @@ def get_hexagonal_prism(edge_length=1., orientation='y',
     l = edge_length
 
     if orientation == 'y':
-        right = XPlane(x0=sqrt(3.)/2.*l)
-        left = XPlane(x0=-sqrt(3.)/2.*l)
+        right = XPlane(x0=sqrt(3.)/2*l, boundary_type=boundary_type)
+        left = XPlane(x0=-sqrt(3.)/2*l, boundary_type=boundary_type)
         c = sqrt(3.)/3.
 
         # y = -x/sqrt(3) + a
@@ -2188,9 +2110,14 @@ def get_hexagonal_prism(edge_length=1., orientation='y',
         prism = -right & +left & -upper_right & -upper_left & \
                 +lower_right & +lower_left
 
+        if boundary_type == 'periodic':
+            right.periodic_surface = left
+            upper_right.periodic_surface = lower_left
+            lower_right.periodic_surface = upper_left
+
     elif orientation == 'x':
-        top = YPlane(y0=sqrt(3.)/2.*l)
-        bottom = YPlane(y0=-sqrt(3.)/2.*l)
+        top = YPlane(y0=sqrt(3.)/2*l, boundary_type=boundary_type)
+        bottom = YPlane(y0=-sqrt(3.)/2*l, boundary_type=boundary_type)
         c = sqrt(3.)
 
         # y = -sqrt(3)*(x - a)
@@ -2207,148 +2134,68 @@ def get_hexagonal_prism(edge_length=1., orientation='y',
         prism = -top & +bottom & -upper_right & +lower_right & \
                             +lower_left & -upper_left
 
+        if boundary_type == 'periodic':
+            top.periodic_surface = bottom
+            upper_right.periodic_surface = lower_left
+            lower_right.periodic_surface = upper_left
+
+    # Handle rounded corners if given
     if corner_radius > 0.:
+        if boundary_type == 'periodic':
+            raise ValueError('Periodic boundary conditions not permitted when '
+                             'rounded corners are used.')
+
+        c = sqrt(3.)/2
+        t = l - corner_radius/c
+
+        # Cylinder with corner radius and boundary type pre-applied
+        cyl1 = partial(ZCylinder, R=corner_radius, boundary_type=boundary_type)
+        cyl2 = partial(ZCylinder, R=corner_radius/(2*c),
+                       boundary_type=boundary_type)
+
         if orientation == 'x':
-            c = sqrt(3.)
-            x_min_y_min_in = ZCylinder(name='x min y min in',
-                                       x0=-0.5  * (l - 2./c*corner_radius),
-                                       y0=-c/2. * (l - 2./c*corner_radius),
-                                       R=corner_radius,
-                                       boundary_type=boundary_type)
-            x_min_y_max_in = ZCylinder(name='x min y max in',
-                                       x0= 0.5  * (l - 2./c*corner_radius),
-                                       y0=-c/2. * (l - 2./c*corner_radius),
-                                       R=corner_radius,
-                                       boundary_type=boundary_type)
-            x_max_y_min_in = ZCylinder(name='x max y min in',
-                                       x0=-0.5  * (l - 2./c*corner_radius),
-                                       y0= c/2. * (l - 2./c*corner_radius),
-                                       R=corner_radius,
-                                       boundary_type=boundary_type)
-            x_max_y_max_in = ZCylinder(name='x max y max in',
-                                       x0= 0.5  * (l - 2./c*corner_radius),
-                                       y0= c/2. * (l - 2./c*corner_radius),
-                                       R=corner_radius,
-                                       boundary_type=boundary_type)
-            x_min_in = ZCylinder(name='x min in',
-                                 x0=-(l - 2./c*corner_radius),
-                                 y0=0.,
-                                 R=corner_radius,
-                                 boundary_type=boundary_type)
-            x_max_in = ZCylinder(name='x max in',
-                                 x0= (l - 2./c*corner_radius),
-                                 y0=0.,
-                                 R=corner_radius,
-                                 boundary_type=boundary_type)
+            x_min_y_min_in = cyl1(name='x min y min in', x0=-t/2, y0=-c*t)
+            x_min_y_max_in = cyl1(name='x min y max in', x0=t/2, y0=-c*t)
+            x_max_y_min_in = cyl1(name='x max y min in', x0=-t/2, y0=c*t)
+            x_max_y_max_in = cyl1(name='x max y max in', x0=t/2, y0=c*t)
+            x_min_in = cyl1(name='x min in', x0=-t)
+            x_max_in = cyl1(name='x max in', x0=t)
 
-            x_min_y_min_out = ZCylinder(name='x min y min out',
-                                        x0=-0.5  * l,
-                                        y0=-c/2. * l,
-                                        R=corner_radius/c,
-                                        boundary_type=boundary_type)
-            x_min_y_max_out = ZCylinder(name='x min y max out',
-                                        x0=+0.5  * l,
-                                        y0=-c/2. * l,
-                                        R=corner_radius/c,
-                                        boundary_type=boundary_type)
-            x_max_y_min_out = ZCylinder(name='x max y min out',
-                                        x0=-0.5  * l,
-                                        y0= c/2. * l,
-                                        R=corner_radius/c,
-                                        boundary_type=boundary_type)
-            x_max_y_max_out = ZCylinder(name='x max y max out',
-                                        x0= 0.5  * l,
-                                        y0= c/2. * l,
-                                        R=corner_radius/c,
-                                        boundary_type=boundary_type)
-            x_min_out = ZCylinder(name='x min out',
-                                  x0=-l,
-                                  y0=0.,
-                                  R=corner_radius/c,
-                                  boundary_type=boundary_type)
-            x_max_out = ZCylinder(name='x max out',
-                                  x0= l,
-                                  y0=0.,
-                                  R=corner_radius/c,
-                                  boundary_type=boundary_type)
+            x_min_y_min_out = cyl2(name='x min y min out', x0=-l/2, y0=-c*l)
+            x_min_y_max_out = cyl2(name='x min y max out', x0=l/2, y0=-c*l)
+            x_max_y_min_out = cyl2(name='x max y min out', x0=-l/2, y0=c*l)
+            x_max_y_max_out = cyl2(name='x max y max out', x0=l/2, y0=c*l)
+            x_min_out = cyl2(name='x min out', x0=-l)
+            x_max_out = cyl2(name='x max out', x0=l)
 
-            corners = (+x_min_y_min_in & -x_min_y_min_out) | \
-                      (+x_min_y_max_in & -x_min_y_max_out) | \
-                      (+x_max_y_min_in & -x_max_y_min_out) | \
-                      (+x_max_y_max_in & -x_max_y_max_out) | \
-                      (+x_min_in       & -x_min_out      ) | \
-                      (+x_max_in       & -x_max_out      )
+            corners = (+x_min_y_min_in & -x_min_y_min_out |
+                       +x_min_y_max_in & -x_min_y_max_out |
+                       +x_max_y_min_in & -x_max_y_min_out |
+                       +x_max_y_max_in & -x_max_y_max_out |
+                       +x_min_in & -x_min_out |
+                       +x_max_in & -x_max_out)
 
         elif orientation == 'y':
-            c = sqrt(3.)
-            x_min_y_min_in = ZCylinder(name='x min y min in',
-                                       x0=-c/2. * (l - 2./c*corner_radius),
-                                       y0=-0.5  * (l - 2./c*corner_radius),
-                                       R=corner_radius,
-                                       boundary_type=boundary_type)
-            x_min_y_max_in = ZCylinder(name='x min y max in',
-                                       x0=-c/2. * (l - 2./c*corner_radius),
-                                       y0= 0.5  * (l - 2./c*corner_radius),
-                                       R=corner_radius,
-                                       boundary_type=boundary_type)
-            x_max_y_min_in = ZCylinder(name='x max y min in',
-                                       x0= c/2. * (l - 2./c*corner_radius),
-                                       y0=-0.5  * (l - 2./c*corner_radius),
-                                       R=corner_radius,
-                                       boundary_type=boundary_type)
-            x_max_y_max_in = ZCylinder(name='x max y max in',
-                                       x0= c/2. * (l - 2./c*corner_radius),
-                                       y0= 0.5  * (l - 2./c*corner_radius),
-                                       R=corner_radius,
-                                       boundary_type=boundary_type)
-            y_min_in = ZCylinder(name='y min in',
-                                 x0=0.,
-                                 y0=-(l - 2./c*corner_radius),
-                                 R=corner_radius,
-                                 boundary_type=boundary_type)
-            y_max_in = ZCylinder(name='y max in',
-                                 x0=0.,
-                                 y0= (l - 2./c*corner_radius),
-                                 R=corner_radius,
-                                 boundary_type=boundary_type)
+            x_min_y_min_in = cyl1(name='x min y min in', x0=-c*t, y0=-t/2)
+            x_min_y_max_in = cyl1(name='x min y max in', x0=-c*t, y0=t/2)
+            x_max_y_min_in = cyl1(name='x max y min in', x0=c*t, y0=-t/2)
+            x_max_y_max_in = cyl1(name='x max y max in', x0=c*t, y0=t/2)
+            y_min_in = cyl1(name='y min in', y0=-t)
+            y_max_in = cyl1(name='y max in', y0=t)
 
-            x_min_y_min_out = ZCylinder(name='x min y min out',
-                                        x0=-c/2. * l,
-                                        y0=-0.5  * l,
-                                        R=corner_radius/c,
-                                        boundary_type=boundary_type)
-            x_min_y_max_out = ZCylinder(name='x min y max out',
-                                        x0=-c/2. * l,
-                                        y0=+0.5  * l,
-                                        R=corner_radius/c,
-                                        boundary_type=boundary_type)
-            x_max_y_min_out = ZCylinder(name='x max y min out',
-                                        x0= c/2. * l,
-                                        y0=-0.5  * l,
-                                        R=corner_radius/c,
-                                        boundary_type=boundary_type)
-            x_max_y_max_out = ZCylinder(name='x max y max out',
-                                        x0= c/2. * l,
-                                        y0= 0.5  * l,
-                                        R=corner_radius/c,
-                                        boundary_type=boundary_type)
-            y_min_out = ZCylinder(name='y min out',
-                                  x0=0.,
-                                  y0=-l,
-                                  R=corner_radius/c,
-                                  boundary_type=boundary_type)
-            y_max_out = ZCylinder(name='y max out',
-                                  x0=0.,
-                                  y0= l,
-                                  R=corner_radius/c,
-                                  boundary_type=boundary_type)
+            x_min_y_min_out = cyl2(name='x min y min out', x0=-c*l, y0=-l/2)
+            x_min_y_max_out = cyl2(name='x min y max out', x0=-c*l, y0=l/2)
+            x_max_y_min_out = cyl2(name='x max y min out', x0=c*l, y0=-l/2)
+            x_max_y_max_out = cyl2(name='x max y max out', x0=c*l, y0=l/2)
+            y_min_out = cyl2(name='y min out', y0=-l)
+            y_max_out = cyl2(name='y max out', y0=l)
 
-            corners = (+x_min_y_min_in & -x_min_y_min_out) | \
-                      (+x_min_y_max_in & -x_min_y_max_out) | \
-                      (+x_max_y_min_in & -x_max_y_min_out) | \
-                      (+x_max_y_max_in & -x_max_y_max_out) | \
-                      (+y_min_in       & -y_min_out      ) | \
-                      (+y_max_in       & -y_max_out      )
+            corners = (+x_min_y_min_in & -x_min_y_min_out |
+                       +x_min_y_max_in & -x_min_y_max_out |
+                       +x_max_y_min_in & -x_max_y_min_out |
+                       +x_max_y_max_in & -x_max_y_max_out |
+                       +y_min_in & -y_min_out |
+                       +y_max_in & -y_max_out)
 
         prism = prism & ~corners
 
diff --git a/src/geometry.F90 b/src/geometry.F90
index 4856159913..5276edb24c 100644
--- a/src/geometry.F90
+++ b/src/geometry.F90
@@ -392,6 +392,7 @@ contains
     real(8) :: v         ! y-component of direction
     real(8) :: w         ! z-component of direction
     real(8) :: norm      ! "norm" of surface normal
+    real(8) :: d         ! distance between point and plane
     real(8) :: xyz(3)    ! Saved global coordinate
     integer :: i_surface ! index in surfaces
     logical :: found     ! particle found in universe?
@@ -532,6 +533,20 @@ contains
         type is (SurfaceZPlane)
           p % coord(1) % xyz(3) = opposite % z0
         end select
+
+      type is (SurfacePlane)
+        select type (opposite => surfaces(surf % i_periodic) % obj)
+        type is (SurfacePlane)
+          ! Get surface normal for opposite plane
+          xyz(:) = opposite % normal(p % coord(1) % xyz)
+
+          ! Determine distance to plane
+          norm = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
+          d = opposite % evaluate(p % coord(1) % xyz) / norm
+
+          ! Move particle along normal vector based on distance
+          p % coord(1) % xyz(:) = p % coord(1) % xyz(:) - d*xyz
+        end select
       end select
 
       ! Reassign particle's surface
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 6c1a40ee3a..b05d560752 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -1658,6 +1658,15 @@ contains
             surf % i_periodic = surface_dict % get_key(surf % i_periodic)
           end if
 
+        type is (SurfacePlane)
+          if (surf % i_periodic == NONE) then
+            call fatal_error("No matching periodic surface specified for &
+                 &periodic boundary condition on surface " // &
+                 trim(to_str(surf % id)) // ".")
+          else
+            surf % i_periodic = surface_dict % get_key(surf % i_periodic)
+          end if
+
         class default
           call fatal_error("Periodic boundary condition applied to &
                &non-planar surface.")