Crystals with complicated geometry are often observed with mixed chemical occupancy among Wyckoff sites, presenting a unique challenge for accurate atomic modeling. Similar systems possessing exact occupancy on all sites can exhibit superstructural ordering, dramatically inflating the unit cell size. Here, we propose a site permutation search (SPS) algorithm for discovering the optimal atomic decoration on fixed crystalline geometries. Our approach relies on the evidence that, for a given chemical composition, a crystal graph convolutional neural network estimates the correct energetic ordering of different atomic decorations, as predicted by electronic structure calculations. This provides a suitable energy landscape, in which site occupation is optimized to discover the chemical decoration of crystals exhibiting mixed or disordered occupancy, or superstructural ordering. The optimization algorithm is based on Monte Carlo moves combined with simulated annealing and basin-hopping techniques. Verification of the procedure is carried out on several known compounds, including the superstructurally ordered clathrate compound Rb8Ga27Sb16 and vacancy-ordered perovskite Cs2SnI6. This strategy provides a fast, accurate, and scalable method for determining favorable decoration in complex crystals.