def get_111_surface(dir, monolayers, vacuum, constrain_layers = 0)

      if dir == "A"
        top_atom = self.anion
        bot_atom = self.cation
      elsif dir == "B"
        top_atom = self.cation
        bot_atom = self.anion
      else
        raise "Direction must be either A or B"
      end
      
      # The atoms on a FCC 
      as1 = Atom.new(0.0, 0.0, 0.0, top_atom)
      ga1 = Atom.new(0.0, 0.0, -sqrt(3)/4*self.lattice_const, bot_atom)

      # The lattice Vectors
      v1 = Vector[0.5*sqrt(2), 0.0, 0.0]*self.lattice_const
      v2 = Vector[sqrt(2)*0.25, sqrt(6)*0.25, 0.0]*self.lattice_const
      v3 = Vector[sqrt(2)*0.25, sqrt(2.0/3.0)*0.25, -1*sqrt(4.0/3.0)*0.5]*self.lattice_const

      # The unit cell
      zb = Geometry.new([as1, ga1], [v1, v2, v3])

      # The Miller Indices
      millerX = [-1, 1, 0]  # Orientation of the crystal pointing in the cartesian +x axis
      millerY = [1, 1, -2]  # Orientation of the crystal pointing in the cartesian +y axis
      millerZ = [-1, -1, -1] # Orientation of the crystal pointing in the cartesian +z axis

      zb.set_miller_indices(millerX, millerY, millerZ)

      # Repeat the unit cell and add vacuum
      if 0 < vacuum 
        # We actually repeat the unit cell monolayers+1 times because
        # I will strip off the top and bottom atoms to make the proper surface
        zb = zb.repeat(1,1,monolayers+1)
        
        bilayerSep = v3[2]
        zb.lattice_vectors[2] = Vector[0, 0, (monolayers-1)*(bilayerSep.abs) + vacuum]

        # Strip off the top and bottom atom
        minZ = zb.atoms.min{|a,b| a.z <=> b.z}.z
        maxZ = zb.atoms.max{|a,b| a.z <=> b.z}.z

        zb.atoms.reject!{|a| a.z == maxZ}
        zb.atoms.reject!{|a| a.z == minZ}

        # Constrain the bottom layers if requested
        if 0 < constrain_layers
          # get the min again because we removed the atoms at minZ above
          minZ = zb.atoms.min{|a,b| a.z <=> b.z}.z
          constrain_below = minZ + bilayerSep.abs*constrain_layers
          zb.atoms.each{|a|
            if (a.z < constrain_below)
              a.constrain = ".true."
            end
          }
        end
      end
      
      zb
    end
def get_111_surface(dir, monolayers, vacuum, constrain_layers = 0)

      if dir == "A"
        top_atom = self.anion
        bot_atom = self.cation
      elsif dir == "B"
        top_atom = self.cation
        bot_atom = self.anion
      else
        raise "Direction must be either A or B"
      end
      
      # The atoms on a FCC 
      as1 = Atom.new(0.0, 0.0, 0.0, top_atom)
      ga1 = Atom.new(0.0, 0.0, -sqrt(3)/4*self.lattice_const, bot_atom)

      # The lattice Vectors
      v1 = Vector[0.5*sqrt(2), 0.0, 0.0]*self.lattice_const
      v2 = Vector[sqrt(2)*0.25, sqrt(6)*0.25, 0.0]*self.lattice_const
      v3 = Vector[sqrt(2)*0.25, sqrt(2.0/3.0)*0.25, -1*sqrt(4.0/3.0)*0.5]*self.lattice_const

      # The unit cell
      zb = Geometry.new([as1, ga1], [v1, v2, v3])

      # The Miller Indices
      millerX = [-1, 1, 0]  # Orientation of the crystal pointing in the cartesian +x axis
      millerY = [1, 1, -2]  # Orientation of the crystal pointing in the cartesian +y axis
      millerZ = [-1, -1, -1] # Orientation of the crystal pointing in the cartesian +z axis

      zb.set_miller_indices(millerX, millerY, millerZ)

      # Repeat the unit cell and add vacuum
      if 0 < vacuum 
        # We actually repeat the unit cell monolayers+1 times because
        # I will strip off the top and bottom atoms to make the proper surface
        zb = zb.repeat(1,1,monolayers+1)
        
        bilayerSep = v3[2]
        zb.lattice_vectors[2] = Vector[0, 0, (monolayers-1)*(bilayerSep.abs) + vacuum]

        # Strip off the top and bottom atom
        minZ = zb.atoms.min{|a,b| a.z <=> b.z}.z
        maxZ = zb.atoms.max{|a,b| a.z <=> b.z}.z

        zb.atoms.reject!{|a| a.z == maxZ}
        zb.atoms.reject!{|a| a.z == minZ}

        # Constrain the bottom layers if requested
        if 0 < constrain_layers
          # get the min again because we removed the atoms at minZ above
          minZ = zb.atoms.min{|a,b| a.z <=> b.z}.z
          constrain_below = minZ + bilayerSep.abs*constrain_layers
          zb.atoms.each{|a|
            if (a.z < constrain_below)
              a.constrain = ".true."
            end
          }
        end
      end
      
      zb
    end
Return a unit cell for a slab of 111
 dir is either "A" or "B" for the cation or anion terminated slab
 specify the number of atomic monolayers
 and the vacuum thickness in angstrom