以前のリビジョンの文書です


波動関数

[kimi@ssp1 ncfiles]$ ./getwfs 
usage: getwfs [-w] input_nc_file output_file_prefix

options:
  -h, --help          show this help message and exit
  -w, --wavefunction  create xyz format files for all wavefunctions

[kimi@ssp1 ncfiles]$ cat ./getdensity ./getwfs 
#!/usr/bin/env python
from sys import exit
from optparse import OptionParser

def check_error(opt, NumberOfBands, NumberOfKpoints, isSpinPolarized):
  error = False
  if (opt.indecies is not None) and (opt.energy is not None):
    print "error: do not specify both options -i and -e"
    error = True
  if opt.indecies is not None:
    (nband, nkpt, nspin) = opt.indecies
    if not (nband < NumberOfBands):
      print "error: nband must be less than %d" % NumberOfBands
      error = True
    if not (nkpt < NumberOfKpoints):
      print "error: nkpt must be less than %d" % NumberOfKpoints 
      error = True
    if isSpinPolarized:
      if not ((nspin == 0) or (nspin == 1)):
        print "error: nspin must be 0 or 1"
        error = True
    else:
      if not (nspin == 0):
        print "error: nspin must be 0"
        error = True
    if error:
      return ('error', nband, nkpt, nspin, 0, 0)
    else:
      return ('index', nband, nkpt, nspin, 0, 0)
  elif (opt.energy is not None):
    (min_energy, max_energy) = opt.energy
    if (min_energy > max_energy):
      (min_energy, max_energy) = (max_energy, min_energy)
    if error:
      return ('error', 0, 0, 0, min_energy, max_energy)
    else:
      return ('energy', 0, 0, 0, min_energy, max_energy)
  else:
    return ('total', 0, 0, 0, 0, 0)


from Dacapo import Dacapo
from ASE.IO.Cube import WriteCube

cmd = OptionParser(usage = '%prog [-i|-e] input_nc_file output_cube_file')
 
cmd.add_option('-i', '--indecies', type = 'int', nargs = 3,
               help = 'density of N-th Band, K-th Kpoint, and S-th Spin state',
               metavar = 'N K S')

cmd.add_option('-e', '--energy', type = 'float', nargs = 2,
               help = 'density of states with energies from E1[eV] to E2[eV]',
               metavar = 'E1 E2')



(opt, argv) = cmd.parse_args()

if len(argv) != 2:
    cmd.print_help()
    raise SystemExit

ncfile = argv[0]
cubefile = argv[1]

model = Dacapo.ReadAtoms(ncfile)
calculator = model.GetCalculator()
NumberOfBands = calculator.GetNumberOfBands()
NumberOfKpoints = len(calculator.GetIBZKPoints())
isSpinPolarized = calculator.GetSpinPolarized()

(commandtype, nband, nkpt, nspin, min_energy, max_energy) = check_error(opt, NumberOfBands, NumberOfKpoints, isSpinPolarized)

if commandtype is 'index':
  wf = calculator.GetWaveFunctionArray(band = nband, kpt = nkpt, spin = nspin)
  density = abs(wf)*abs(wf)
  WriteCube(model, density, cubefile)
elif commandtype is 'energy':
  if (isSpinPolarized):
    NumberOfSpins = 2
  else:
    NumberOfSpins = 1
  density = calculator.GetDensityArray()*0.0
  out_of_range = True
  min_range = 0.0
  max_range = 0.0
  for nspin in range(NumberOfSpins):
    for nkpt in range(NumberOfKpoints):
      energies = calculator.GetEigenvalues(kpt = nkpt, spin = nspin)
      for nband in range(NumberOfBands):
        if (energies[nband] < min_range):
          min_range = energies[nband] 
        if (energies[nband] > max_range):
          max_range = energies[nband]
        if (min_energy <= energies[nband]) and (energies[nband] <= max_energy):
          out_of_range = False
          wf = calculator.GetWaveFunctionArray(band = nband, kpt = nkpt, spin = nspin)
          density = density + abs(wf)*abs(wf)
  if out_of_range:
    print "error: energy range should be from %f to %f" % (min_range, max_range)
    exit()
  WriteCube(model, density, cubefile)
else:
  density = calculator.GetDensityArray()
  WriteCube(model, density, cubefile)

#!/usr/bin/env python
from math import log10
from optparse import OptionParser

from Dacapo import Dacapo
from ASE.IO.Cube import WriteCube

cmd = OptionParser(usage = '%prog [-w] input_nc_file output_file_prefix')
cmd.add_option('-w', '--wavefunction', action = "store_true", default = False,
               help = 'create xyz format files for all wavefunctions')


(opt, argv) = cmd.parse_args()

if len(argv) != 2:
    cmd.print_help()
    raise SystemExit

ncfile = argv[0]
cubefileprefix = argv[1]

model = Dacapo.ReadAtoms(ncfile)
calculator = model.GetCalculator()

if opt.wavefunction is not True:
  density = calculator.GetDensityArray()
  cubefile = '%s.xyz' % cubefileprefix
  WriteCube(model, density, cubefile)
else:
  NumberOfBands = calculator.GetNumberOfBands()
  NumberOfKpoints = len(calculator.GetIBZKPoints())
  isSpinPolarized = calculator.GetSpinPolarized()
  if (isSpinPolarized):
    NumberOfSpins = 2
  else:
    NumberOfSpins = 1

  for nband in range(NumberOfBands):
    for nkpt in range(NumberOfKpoints):
      for nspin in range(NumberOfSpins):
        cubefile = '%s' % cubefileprefix
        fmts = '%%0%dd' % (log10(NumberOfBands) + 1)
        cubefile += fmts % nband
        if (NumberOfKpoints > 1):
          fmts = 'k%%0%dd' % (log10(NumberOfKpoints) + 1)
          cubefile += fmts % nkpt
        if (NumberOfSpins > 1):
          cubefile += 's%1d' % nspin
        cubefile += '.cube'
        print cubefile
        wf = calculator.GetWaveFunctionArray(band = nband,
                                              kpt = nkpt,
                                             spin = nspin)
        WriteCube(model, wf, cubefile)


[kimi@ssp1 ncfiles]$ 
ab_initio/波動関数.1233416997.txt.gz · 最終更新: 2009/02/01 00:49 by kimi
www.chimeric.de Creative Commons License Valid CSS Driven by DokuWiki do yourself a favour and use a real browser - get firefox!! Recent changes RSS feed Valid XHTML 1.0