def get_001_surface(monolayers, vacuum, constrain_layers = 0)
      anion = Atom.new(0,0,0,self.cation)
      cation = Atom.new(0.25*self.lattice_const, 0.25*self.lattice_const, 0.25*self.lattice_const, self.anion)
      v1 = Vector[0.5, 0.5, 0]*self.lattice_const
      v2 = Vector[-0.5,0.5,0]*self.lattice_const
      v3 = Vector[0.5, 0, 0.5]*self.lattice_const
      
      zb = Geometry.new([anion, cation], [v1,v2,v3])
      millerX = [1,0,0]
      millerY = [0,1,0]
      millerZ = [0,0,1]
      zb.set_miller_indices(millerX, millerY, millerZ)
      
      # Repeat the unit cell.  The unit cell is a bi-layer so divide by 2
      zb = zb.repeat(1,1,(monolayers/2).ceil)

      if 0 < vacuum
        # Add vacuum
        monolayerSep = v3[2]/2
        zb.lattice_vectors[2] = Vector[0, 0, (monolayers-1)*monolayerSep.abs + vacuum.to_f]
        # Move everything into a nice tidy unit cell. 
        zb = zb.correct
      end

      minZ = zb.atoms.min{|a,b| a.z <=> b.z}.z
      
      # Reject the top layer of atoms if an odd number of monolayers was requested.
      # This is necessary because the primitive cell is a bilayer
      zb.atoms.reject! {|a| 
        a.z >= (minZ + monolayerSep.abs*monolayers)
      }
      
      # Constrain the bottom layers
      zb.atoms.each{|a|
        if (a.z < minZ + monolayerSep.abs*constrain_layers)
          a.constrain = ".true."
        end
      }
      
      # Return the completed unit cell
      return zb
    end
def get_001_surface(monolayers, vacuum, constrain_layers = 0)
      anion = Atom.new(0,0,0,self.cation)
      cation = Atom.new(0.25*self.lattice_const, 0.25*self.lattice_const, 0.25*self.lattice_const, self.anion)
      v1 = Vector[0.5, 0.5, 0]*self.lattice_const
      v2 = Vector[-0.5,0.5,0]*self.lattice_const
      v3 = Vector[0.5, 0, 0.5]*self.lattice_const
      
      zb = Geometry.new([anion, cation], [v1,v2,v3])
      millerX = [1,0,0]
      millerY = [0,1,0]
      millerZ = [0,0,1]
      zb.set_miller_indices(millerX, millerY, millerZ)
      
      # Repeat the unit cell.  The unit cell is a bi-layer so divide by 2
      zb = zb.repeat(1,1,(monolayers/2).ceil)

      if 0 < vacuum
        # Add vacuum
        monolayerSep = v3[2]/2
        zb.lattice_vectors[2] = Vector[0, 0, (monolayers-1)*monolayerSep.abs + vacuum.to_f]
        # Move everything into a nice tidy unit cell. 
        zb = zb.correct
      end

      minZ = zb.atoms.min{|a,b| a.z <=> b.z}.z
      
      # Reject the top layer of atoms if an odd number of monolayers was requested.
      # This is necessary because the primitive cell is a bilayer
      zb.atoms.reject! {|a| 
        a.z >= (minZ + monolayerSep.abs*monolayers)
      }
      
      # Constrain the bottom layers
      zb.atoms.each{|a|
        if (a.z < minZ + monolayerSep.abs*constrain_layers)
          a.constrain = ".true."
        end
      }
      
      # Return the completed unit cell
      return zb
    end
Return a unit cell for a slab of 001
 Specify the number of atomic monolayers,
 the vacuum thickness in angstrom,
 and the number of layers to constrain at the base of the slab