The multipole boundary conditions is implemented using a real space multigrid poisson solver. In this case there is no periodic boundary conditions, we are simulating a truly isolated molecule. This is a standard approach in molecular quantum chemistry softwares, but maybe a less know methodology in the solid state community.
The potential at the boundary is determined by calculating the monopole, dipole and quadrupole moments of the charge distribution inside the box, and using these moments to extrapolate the value of the electro-static potential at the boundary of the box. Thus, the asymptotic values are determined using the charge multipoles, and we only need to describe the potential inside the box, using the asymptotic values as boundary conditions.
Currently, the main approximation is that we use a first order stensil for the multigrid solver, we are working on also implementing higher order stensils.
The reference you give is not relevant.
To get the box size, I recommend that you calculate the ElectrostaticDifferencePotential, and get the box using
potential = ElectrostaticDifferencePotential(molecule_configuration)
a, b, c = potential.unitCell()
print a, b, c