1 from Dacapo import Dacapo
   2 from ASE import Atom
   3 from ASE.IO.Cube import WriteCube
   4 inputfile='yourfile.nc'
   5 atoms0=Dacapo.ReadAtoms(inputfile)
   6 calc0=atoms0.GetCalculator()
   7 sys_density=calc0.GetDensityArray()
   8 
   9 atoms1=Dacapo.ReadAtoms(inputfile)
  10 #remove the slab
  11 del atoms1[0:2]
  12 calc1 = Dacapo(nbands=2 * 5,
  13               kpts=(4, 4, 1),
  14               planewavecutoff=340,
  15               densitycutoff=400)
  16 calc1.SetTxtFile('atom.txt')
  17 calc1.SetNetCDFFile('atom.nc')
  18 
  19 atoms1.SetCalculator(calc1)
  20 atoms1.GetPotentialEnergy()
  21 atom_density=calc1.GetDensityArray()
  22 
  23 atoms2=Dacapo.ReadAtoms(inputfile)
  24 #remove the ad atom
  25 del atoms2[2]
  26 calc2 = Dacapo(nbands=2 * 5,
  27               kpts=(4, 4, 1),
  28               planewavecutoff=340,
  29               densitycutoff=400)
  30 calc2.SetTxtFile('slab.txt')
  31 calc2.SetNetCDFFile('slab.nc')
  32 
  33 atoms2.SetCalculator(calc2)
  34 atoms2.GetPotentialEnergy()
  35 slab_density=calc2.GetDensityArray()
  36 
  37 density_diff=(sys_density-atom_density-slab_density)
  38 
  39 
  40 WriteCube(atoms0,density_diff,'densitydiff_array.cube')