def METHOD_NAME(els, density, sx=None, sy=None):
    if isinstance(els, str):
        syms = string2symbols(els)
    else:
        syms = ""
        for sym, n in els:
            syms += n*sym
    r = np.random.rand(len(syms), 3)
    a = ase.Atoms(syms, positions=r, cell=[1,1,1], pbc=True)
    mass = np.sum(a.get_masses())
    if sx is not None and sy is not None:
        sz = ( 1e24*mass/(density*mol) )/(sx*sy)
        a.set_cell([sx,sy,sz], scale_atoms=True)
    else:
        a0 = ( 1e24*mass/(density*mol) )**(1./3)
        a.set_cell([a0,a0,a0], scale_atoms=True)
    a.set_initial_charges([0]*len(a))
    return a