1 from Dacapo import Dacapo
   2 from ASE import Atom, ListOfAtoms
   3 from ASE.Visualization.VMD import VMD
   4 from ASE.Visualization.gnuplot import gnuplot
   5 
   6 bulk = ListOfAtoms([Atom('Al', (0, 0, 0))] )
   7 a = 4.05
   8 b = a / 2
   9 bulk.SetUnitCell([(0, b, b),
  10                   (b, 0, b),
  11                   (b, b, 0)])
  12 
  13 bulk_plot = VMD(bulk, repeat=(5,5,5))
  14 
  15 calc = Dacapo(kpts=(1, 1, 1),          # set the k-points (Monkhorst-Pack)
  16               planewavecutoff=200,     # planewavecutoff in eV
  17               nbands=6,                # set the number of electronic bands
  18               usesymm=True,            # use symmetry to reduce the k-point set
  19               out='Al-fcc.nc',         # define the out netcdf file
  20               txtout='Al-fcc.txt')     # define the ascii out file
  21 
  22 bulk.SetCalculator(calc)
  23 calc.SetDensityCutoff(400)
  24 
  25 # make a plot of the convergence with respect to k-points
  26 kpt_energies = []
  27 for n in [1, 2, 4, 6]: 
  28     calc.SetBZKPoints((n, n, n))
  29     energy = bulk.GetPotentialEnergy() 
  30     kpt_energies.append((n, energy))
  31 
  32 kpt_plot = gnuplot(kpt_energies) 
  33 
  34 
  35 # make a plot of the convergence with respect to planewavecutoff
  36 n = 2
  37 calc.SetBZKPoints((n, n, n))
  38 pw_energies = []
  39 for planewavecutoff in [75, 100, 150, 200, 250, 300]:
  40     calc.SetPlaneWaveCutoff(planewavecutoff)
  41     energy = bulk.GetPotentialEnergy()
  42     pw_energies.append((planewavecutoff, energy))
  43 
  44 pw_plot = gnuplot(pw_energies)