1 #!/usr/bin/env python
   2 
   3 from math import sqrt
   4 
   5 import Numeric as num
   6 
   7 from Dacapo import Dacapo
   8 from ASE.Filters.Subset import Subset
   9 from ASE.Filters.NudgedElasticBand import NudgedElasticBand
  10 from ASE.Dynamics.QuasiNewton import QuasiNewton
  11 from ASE.Trajectories.NetCDFTrajectory import NetCDFTrajectory
  12 from ASE import Atom,ListOfAtoms
  13 
  14 # setup the static Al slab 
  15 alslab =  ListOfAtoms([ Atom('Al',(0, 0, 0)),
  16                         Atom('Al',(0, 0.5, 0)),
  17                         Atom('Al',(0.5, 0, 0)),
  18                         Atom('Al',(0.5, 0.5, 0)),
  19                         Atom('Al',(0.3333333, 0.166666667, 0.125*2.)),
  20                         Atom('Al',(0.3333333, 0.666666667, 0.125*2.)),
  21                         Atom('Al',(0.8333333, 0.166666667, 0.125*2.)),
  22                         Atom('Al',(0.8333333, 0.666666667, 0.125*2.))])
  23 
  24 initial = alslab.Copy()
  25 initial.append(Atom('O',(0.166586,0.33341, -0.043448*2.),tag=1))
  26 
  27 final = alslab.Copy()
  28 final.append(Atom('O',(0.33331, 0.166677, -0.04417119*2. ),tag=1))
  29 
  30 unitcell = [[5.72756492761103, 0, 0], 
  31             [-2.86378246380552, 4.96021672913593, 0], 
  32             [0,0,18.7061487217439/2.]]
  33 
  34 initial.SetUnitCell(unitcell) 
  35 final.SetUnitCell(unitcell)  
  36 
  37 mask=[a.GetTag()==1 for a in initial]
  38 
  39 calc = Dacapo(planewavecutoff = 340, 
  40               densitycutoff = 500,
  41               nbands = 22) 
  42 
  43 calc.StayAliveOff()
  44 calc.SetNetCDFFile('initial.nc') 
  45 initial.SetCalculator(calc) 
  46 initialenergy = initial.GetPotentialEnergy() 
  47 # Create a qn object:
  48 subset = Subset(initial,mask) 
  49 relax = QuasiNewton(subset,fmax=0.05)
  50 relax.Converge(maxsteps=1)
  51 initialenergy1 = initial.GetPotentialEnergy() 
  52 
  53 calc.SetNetCDFFile('final.nc') 
  54 final.SetCalculator(calc) 
  55 finalenergy = final.GetPotentialEnergy() 
  56 # Create a qn object:
  57 subset = Subset(final,mask) 
  58 relax = QuasiNewton(subset,fmax=0.05)
  59 relax.Converge(maxsteps=1)
  60 finalenergy1 = final.GetPotentialEnergy() 
  61 
  62 print 'initial state energy             : ',initialenergy
  63 print 'initial state energy (minimized) : ',initialenergy1
  64 print 'final state energy               : ',finalenergy
  65 print 'final state energy (minimized)   : ',finalenergy1
  66 
  67 atomslist = [Dacapo.ReadAtoms('initial.nc')]
  68 configs   = [Subset(atomslist[0], mask=mask)]
  69 for n in range(2): 
  70     atoms = Dacapo.ReadAtoms('initial.nc') 
  71     atomslist.append(atoms)
  72     configs.append(Subset(atoms, mask=mask))
  73 
  74 atoms = Dacapo.ReadAtoms('final.nc') 
  75 configs.append(Subset(atoms, mask=mask))
  76 atomslist.append(atoms)
  77 
  78 # setup each calculator
  79 for n in range(len(atomslist)): 
  80     calc = atomslist[n].GetCalculator() 
  81     calc.StayAliveOff()
  82     calc.SetTxtFile('out.'+str(n)+'.txt') 
  83     calc.SetNetCDFFile('out.'+str(n)+'.nc') 
  84 
  85 band = NudgedElasticBand(configs)
  86 band.SetInterpolatedPositions()
  87 
  88 # Create a qn object:
  89 relax = QuasiNewton(band,fmax=0.05,forcemin=False)
  90 
  91 # Create a trajectory for the each image
  92 listoftraj = []
  93 for n in range(len(atomslist)): 
  94     path = NetCDFTrajectory('image'+str(n)+'.nc', atomslist[n])
  95     print 'n= ',n,atomslist[n].GetPotentialEnergy()
  96     listoftraj.append(path) 
  97     path.Update()
  98     relax.Attach(path) 
  99 
 100 
 101 relax.Converge()
 102 
 103 # make a trajectory of the final configurations in neb-path.nc
 104 atomslist = [listoftraj[n].GetListOfAtoms(-1) for n in range(len(listoftraj))]
 105 
 106 atom0 = atomslist[0]
 107 newtraj = NetCDFTrajectory('neb-path.nc',atom0) 
 108 newtraj.Update()
 109 
 110 for atom in atomslist[1:]: 
 111     atom0.SetCartesianPositions(atom.GetCartesianPositions())
 112     atom0.SetCalculator(atom.GetCalculator())
 113     newtraj.Update()