====== gettotaldensity ======

<code python>
#!/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)
</code>
===== Usage =====
<code sh>
$ gettotaldensity 
usage: gettotaldensity input_nc_file output_cube_file

options:
  -h, --help  show this help message and exit
$ 
</code>