1 from ASE import ListOfAtoms, Atom
   2 from ASE.Dynamics.MDMin import MDMin
   3 from ASE.Dynamics.ConjugateGradient import ConjugateGradient
   4 from ASE.Trajectories import WriteNetCDFTrajectory
   5 from ASE.Filters.Subset import Subset
   6 from ASE.Filters.NudgedElasticBand import NudgedElasticBand
   7 from math import sqrt
   8 
   9 import Asap
  10 
  11 
  12 a = 4.0614
  13 b = a / sqrt(2)
  14 h = b / 2
  15 initial = ListOfAtoms([Atom('Al', (0, 0, 0)),
  16                        Atom('Al', (a / 2, b / 2, -h))],
  17                       periodic=(1, 1, 0),
  18                       cell=(a, b, -2 * h))
  19 initial = initial.Repeat((4, 5, 2))
  20 initial.append(Atom('Al', (a / 2, b / 2, h)))
  21 initial = Asap.ListOfAtoms(initial)
  22 
  23 final = initial.Copy()
  24 final[-1].SetCartesianPosition((a / 2, 3 * b / 2, h))
  25 
  26 # Construct a list of images:
  27 images = [initial]
  28 for i in range(??):
  29     images.append(initial.Copy())
  30 images.append(final)
  31 
  32 # Let all images use an Asap-EMT calculator:
  33 for image in images:
  34     image.SetCalculator(Asap.EMT())
  35 
  36 # Make a mask of zeros and ones that select the dynamic atoms (the
  37 # three topmost layers):
  38 mask = [atom.GetCartesianPosition()[2] > -2.5 * h for atom in initial]
  39 
  40 # The subsets contain only the three topmost layers - the bottom layer
  41 # is filtered out (fixed):
  42 image_subsets = [Subset(image, mask=mask) for image in images]
  43 
  44 # Relax the initial state:
  45 cg = ConjugateGradient(image_subsets[0], fmax=0.02)
  46 cg.Converge()
  47     
  48 # Relax the final state:
  49 cg = ConjugateGradient(image_subsets[-1], fmax=0.02)
  50 cg.Converge()
  51 
  52 # Create a Nudged Elastic Band:
  53 neb = NudgedElasticBand(image_subsets)
  54 
  55 # Mak a starting guess for the minimum energy path (a straight line
  56 # from the initial to the final state):
  57 neb.SetInterpolatedPositions()
  58 
  59 # Use MDMin to relax the path:
  60 minimizer = MDMin(neb, dt=0.08)
  61 
  62 # Relax the path and print out the energy barrier:
  63 for i in range(5):
  64     minimizer.Run(10)
  65     energies = neb.GetListOfEnergies()
  66     Ea = max(energies) - energies[0]
  67     print i, Ea
  68 
  69 # Write the path to a trajectory:
  70 WriteNetCDFTrajectory(images, 'jump1.traj')