gettotaldensity

#!/usr/bin/env python
import string
from math import log10
from optparse import OptionParser
 
from Dacapo import Dacapo
from ASE.IO.Cube import WriteCube
 
cmd = OptionParser(usage = '%prog input_nc_file output_cube_file')
 
(opt, argv) = cmd.parse_args()
 
if len(argv) != 2:
    cmd.print_help()
    raise SystemExit
 
ncfile = argv[0]
cubefile = argv[1]
basename = string.split(cubefile, '.cube')
 
model = Dacapo.ReadAtoms(ncfile)
calculator = model.GetCalculator()
 
if not calculator.GetSpinPolarized():
  filename_total= basename[0] + '.cube'
else:
  filename_total= basename[0] + '_total.cube'
  filename_diff= basename[0] + '_diff.cube'
  filename_spin0= basename[0] + '_spin0.cube'
  filename_spin1= basename[0] + '_spin1.cube'
  density0 = calculator.GetDensityArray(spin = 0)
  density1 = calculator.GetDensityArray(spin = 1)
  WriteCube(model, density0, filename_spin0)
  WriteCube(model, density1, filename_spin1)
  diff = density1 - density0
  WriteCube(model, diff, filename_diff)
 
density = calculator.GetDensityArray()
WriteCube(model, density, filename_total)

$ gettotaldensity 
usage: gettotaldensity input_nc_file output_cube_file
 
options:
  -h, --help  show this help message and exit
$