  We present an approach to identify topological order based on unbiased
infinite projected entangled-pair states (iPEPS) simulations, i.e. where we do
not impose a virtual symmetry on the tensors during the optimization of the
tensor network ansatz. As an example we consider the ground state of the toric
code model in a magnetic field exhibiting $Z_2$ topological order. The
optimization is done by an efficient energy minimization approach based on a
summation of tensor environments to compute the gradient. We show that the
optimized tensors, when brought into the right gauge, are approximately $Z_2$
symmetric, and they can be fully symmetrized a posteriori to generate a stable
topologically ordered state, yielding the correct topological entanglement
entropy and modular S and U matrices. To compute the latter we develop a
variant of the corner-transfer matrix method which is computationally more
efficient than previous approaches based on the tensor renormalization group.
