1 #!/usr/bin/env python
   2 
   3 from Dacapo import Dacapo
   4 from ASE.Visualization.VTK import VTKPlotArray, VTKPlotAtoms
   5 
   6 atoms_dimer = Dacapo.ReadAtoms("NMA-dimer.nc")
   7 atoms_molecule1 = Dacapo.ReadAtoms("NMA_freeze1.nc")
   8 atoms_molecule2 = Dacapo.ReadAtoms("NMA_freeze2.nc")
   9 
  10 calc_dimer = atoms_dimer.GetCalculator()
  11 calc_molecule1 = atoms_molecule1.GetCalculator()
  12 calc_molecule2 = atoms_molecule2.GetCalculator()
  13 
  14 
  15 dens_dimer = calc_dimer.GetDensityArray()
  16 dens_molecule1 = calc_molecule1.GetDensityArray()
  17 dens_molecule2 = calc_molecule2.GetDensityArray()
  18 
  19 dens_diff = dens_dimer - (dens_molecule1 + dens_molecule2)
  20 
  21 plot = VTKPlotArray(dens_diff, atoms_dimer.GetUnitCell())
  22 plot.SetTranslation((-15,-50,0))
  23 
  24 atomplot = VTKPlotAtoms(atoms_dimer, parent = plot)
  25 atomplot.RemoveAvatar(atomplot.unitcell)
  26 plot.Render()
  27 
  28 
  29 #If one wants the atoms white
  30 #atomavatar.SetAtomColors((1,1,1))
  31 
  32 #Get a handle on the isosurfaces and only look at two
  33 #at +-15
  34 isosurface1=plot.GetIsoSurface1()
  35 isosurface1.SetContourRange((-15.1,15))
  36 isosurface1.SetColorRange((-15.1,15))
  37 isosurface1.SetNumberOfContours(2)
  38 plot.Render()
  39 
  40 #The isosurface at -15 is colored blue, the one at +15 yellow
  41 plot.GetColorTable().SetRGBColors(['blue','yellow'])
  42 plot.GetAvatars()[0].GetActorProperty().SetOpacity(0.35)
  43 plot.Render()