Constructing potential energy surface of a periodic system, how to deal with atoms moving out of the supercell?
I think that it is best to consider the calculation of the potential energy (respecting the periodic boundary conditions) as a separate matter from keeping the atoms within the supercell. For simplicity I'm going to assume a cuboidal box with sides $(L_x,L_y,L_z)$, with periodicity in all three Cartesian directions.

Some people like to put the atoms back into the primary cell whenever they move out. This is done by replacing the atomic position vector $(x,y,z)$ by $(x,y,z)+(n_xL_x,n_yL_y,n_zL_z)$, where $(n_x,n_y,n_z)$ are integers, chosen to make the new coordinate values lie in the desired range, often $(\pm\frac{1}{2}L_x,\pm\frac{1}{2}L_y,\pm\frac{1}{2}L_z)$. If you do this every MD step, the integers will usually just need to take the values $0,\pm1$. Other people prefer not to apply this correction, and just keep track of the copy of an atom that diffuses into the periodic-image supercells.

Even if you do this correction to the atom positions, you still need to apply a similar correction every time you calculate a vector between atoms, for the purposes of computing two-body, three-body, potentials and forces. So, a properly written calculator for these quantities will not rely on the atoms lying within the same supercell. If the atomic coordinates are allowed to evolve freely, without periodic boundary corrections, then the integers $(n_x,n_y,n_z)$ will take whatever values are needed to bring the corrected interatomic vector into "minimum image" form, with components in the range $(\pm\frac{1}{2}L_x,\pm\frac{1}{2}L_y,\pm\frac{1}{2}L_z)$. There are standard prescriptions for doing this.

Often, the energy/force calculation is speeded up by a neighbour list calculation, which involves assigning each atom to a subcell of the original supercell. On the face of it, this seems to require that we put the atoms back into the supercell, before doing the assignment. But actually, a properly written function to identify nearest neighbours will have this mapping built in to it, and may just work with copies of the original coordinates. So, again, there should be no need to put the atoms back into the original supercell "by hand", if the neighbour identification is done this way, and if the vectors between neighbours are corrected as just described.

My impression, from a cursory look at ase, is that this mapping is indeed built in to the functions used to identify neighbours "within range". There is a simple Lennard-Jones calculator, for instance, which seems to rely on the neighbour list function for this purpose. The relevant lines are 

        neighbors, offsets = self.nl.get_neighbors(a1)
        cells = np.dot(offsets, cell)
        d = positions[neighbors] + cells - positions[a1]


where I think that offsets contains the integers $(n_x,n_y,n_z)$ and I think that cell contains the cell lengths $(L_x,L_y,L_z)$ in the form of a $3\times3$ matrix (just the diagonal elements would be nonzero in this case).

If you are not using such a function to locate neighbours "within range" and pre-calculate the necessary minimum image corrections to the interatomic vectors, you will need to put it into your calculation of pair and three-body energies and forces, somewhere, yourself. As I mentioned, this is not the same as putting the atoms back into the supercell, when they move out.

Caveat! I am not an expert in ase, and if you want specific advice on that package, you should look at the package documentation, and at simple examples of calculators provided within it. Generally speaking, I believe that advice on debugging programs, and on specific packages, is considered off-topic for this site, so I think that's as far as I should go in any case.