====== makevestafile ======

<code python>
#!/usr/bin/env python
from optparse import OptionParser
from math import sqrt, acos, pi
 
from Dacapo import Dacapo
from ASE.ChemicalElements import Element

def getabs(v):
  _av = sqrt(v[0]**2 + v[1]**2 + v[2]**2)
  return _av

def getangle(v1, v2):
  _av1 = getabs(v1)
  _av2 = getabs(v2)
  _ang = acos((v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2])/(_av1*_av2))*180/pi
  return _ang

def WriteVESTA(filename, atoms, cubefile = None):
  natoms = len(atoms)
  (a1, a2, a3) = atoms.GetUnitCell()
  f = open(filename, 'w')
  f.write("#VESTA_FORMAT_VERSION 2\n")
  f.write('\n')
  f.write("TITLE\n")
  f.write('\n')
  if cubefile is not None:
    f.write("IMPORT_DENSITY\n")
    f.write("+1.000000 %s\n" % cubefile)
  f.write('\n')
  f.write("GROUP\n")
  f.write("1 1 P 1  \n")
  f.write("SYMOP\n")
  f.write(" 0.000000  0.000000  0.000000   1  0  0   0  1  0   0  0  1\n")
  f.write(" -1.0 -1.0 -1.0  0 0 0  0 0 0  0 0 0\n")
  f.write("TRANM\n")
  f.write(" 0.000000  0.000000  0.000000  1  0  0   0  1  0   0  0  1\n")
  f.write("CELLP\n")
  f.write("%10.6f "*6 % (getabs(a1), getabs(a2), getabs(a3), getangle(a2, a3), getangle(a3, a1), getangle(a1, a2)))
  f.write('\n')
  f.write("  0.000000   0.000000   0.000000   0.000000   0.000000   0.000000 \n")
  f.write("STRUC\n")
  for i in range(natoms):
    (x, y, z) = atoms.GetScaledPositions()[i]
    symb = Element(atoms.GetAtomicNumbers()[i]).symbol
    f.write("%3d %2s       H  1.0000 %10.6f %10.6f %10.6f    1a       1\n" % (i + 1, symb, x, y, z))
    f.write("                         0.000000   0.000000   0.000000  0.00\n")
  f.write("  0 0 0 0 0 0 0\n")


cmd = OptionParser(usage = '%prog [-e embeded_cube_file] input_nc_file output_vesta_file')
cmd.add_option('-e', '--embed', type = 'string', nargs = 1,
               help = 'Embed the embeded_cube_file in the output_vesta_file.',
               metavar = 'embeded_cube_file')
 
(opt, argv) = cmd.parse_args()
 
if len(argv) != 2:
    cmd.print_help()
    raise SystemExit
 
ncfile = argv[0]
vestafile = argv[1]
 
model = Dacapo.ReadAtoms(ncfile)
WriteVESTA(vestafile, atoms = model, cubefile = opt.embed)
</code>