from Dacapo import Dacapo
from ASE import Atom
from ASE.IO.Cube import WriteCube
inputfile='yourfile.nc'
atoms0=Dacapo.ReadAtoms(inputfile)
calc0=atoms0.GetCalculator()
sys_density=calc0.GetDensityArray()

atoms1=Dacapo.ReadAtoms(inputfile)
#remove the slab
del atoms1[0:2]
calc1 = Dacapo(nbands=2 * 5,
              kpts=(4, 4, 1),
              planewavecutoff=340,
              densitycutoff=400)
calc1.SetTxtFile('atom.txt')
calc1.SetNetCDFFile('atom.nc')

atoms1.SetCalculator(calc1)
atoms1.GetPotentialEnergy()
atom_density=calc1.GetDensityArray()

atoms2=Dacapo.ReadAtoms(inputfile)
#remove the ad atom
del atoms2[2]
calc2 = Dacapo(nbands=2 * 5,
              kpts=(4, 4, 1),
              planewavecutoff=340,
              densitycutoff=400)
calc2.SetTxtFile('slab.txt')
calc2.SetNetCDFFile('slab.nc')

atoms2.SetCalculator(calc2)
atoms2.GetPotentialEnergy()
slab_density=calc2.GetDensityArray()

density_diff=(sys_density-atom_density-slab_density)


WriteCube(atoms0,density_diff,'densitydiff_array.cube')