1 """ Wannier orbitals for the ethylene molecule 
   2 """
   3 
   4 from ASE import Atom, ListOfAtoms
   5 from Dacapo import Dacapo
   6 from ASE.Utilities.Wannier.Wannier import Wannier
   7 import os
   8 
   9 
  10 # check if we have already calculated the nc file
  11 if not os.path.isfile('ethylene.nc'):
  12 
  13     # ethylene molecule
  14     a = 6.0  # Size of unit cell (Angstrom)
  15     atoms = ListOfAtoms([
  16                         Atom('H', (-1.235, 0.936 , 0 ),tag=0),
  17                         Atom('H', ( 1.235,-0.936 , 0 ),tag=1),
  18                         Atom('H', ( 1.235, 0.936 , 0 ),tag=1),
  19                         Atom('H', (-1.235,-0.936 , 0 ),tag=1),
  20                         Atom('C', ( 0.660, 0.000 , 0 ),tag=1),
  21                         Atom('C', (-0.660, 0.000 , 0 ),tag=1)],
  22                         cell=(a,a,a), periodic=True)
  23 
  24 
  25     pos = atoms.GetCartesianPositions() + 3
  26     atoms.SetCartesianPositions(pos)
  27 
  28     # Dacapo calculator:
  29     calc = Dacapo(planewavecutoff=400, nbands=8, xc='PBE',out='ethylene.nc')
  30     atoms.SetCalculator(calc)
  31     tot = atoms.GetPotentialEnergy()
  32 
  33 atoms = Dacapo.ReadAtoms('ethylene.nc')
  34 calc = atoms.GetCalculator()
  35 
  36 # Initialize the Wannier class
  37 wannier = Wannier(numberofwannier=6,calculator=calc) 
  38 
  39 # Perform the localization with specified convergence criterion
  40 wannier.Localize(tolerance=1.0e-8)
  41 
  42 # Read the centers and radii of the WFs
  43 for center in  wannier.GetCenters(): 
  44     print center
  45 
  46 centers = wannier.GetCentersAsAtoms()
  47 
  48 # plot centers together with atoms
  49 try:
  50     # first try with VMD
  51     from ASE.Visualization.VMD import VMD
  52     VMD(atoms,centers)
  53 except:
  54     from ASE.Visualization.RasMol import RasMol
  55     centers.extend(atoms)
  56     rasmol = RasMol(centers)
  57     
  58 
  59 # make a 3D isosurface plot using VTK
  60 from ASE.Visualization.VTK import VTKPlotElectronicState,VTKPlotAtoms
  61 state = wannier.GetElectronicState(0)
  62 
  63 plot = VTKPlotElectronicState(state)
  64 plot.SetRepresentationToIsoSurface2([0.1])
  65 atomsplot = VTKPlotAtoms(atoms,parent=plot) 
  66 plot.Update()
  67 atomsplot.RemoveAvatar(atomsplot.GetAvatars()[0])
  68 plot.Render()