1 #!/usr/bin/env python
   2 
   3 """ restart a NEB calculations 
   4 
   5 Script to restart a NEB calculation. 
   6 Assume that the configurations are in out.<config-no>.nc, 
   7 config-no = 0,M-1
   8 """ 
   9 
  10 from math import sqrt
  11 
  12 import Numeric as num
  13 
  14 from Dacapo import Dacapo
  15 from ASE.Filters.Subset import Subset
  16 from ASE.Filters.NudgedElasticBand import NudgedElasticBand
  17 from ASE.Dynamics.QuasiNewton import QuasiNewton
  18 from ASE.Trajectories.NetCDFTrajectory import NetCDFTrajectory
  19 from ASE import Atom,ListOfAtoms
  20 
  21 # number of frames in total including end points
  22 M = 4
  23 
  24 # create file names 
  25 filenames = ['initial.nc']
  26 filenames.extend(['out.'+str(n)+'.nc' for n in range(1,M-1)])
  27 filenames.append('final.nc')
  28 
  29 # atoms list 
  30 atomslist = [Dacapo.ReadAtoms(file,save_memory=True) for file in filenames]
  31 
  32 # setup each calculator
  33 # dacapo does not remember StayAliveOff and filenames
  34 for n in range(len(atomslist)): 
  35     calc = atomslist[n].GetCalculator() 
  36     # calc.StayAliveOff()
  37     calc.SetTxtFile('out-restart.'+str(n)+'.txt') 
  38     calc.SetNetCDFFile('out-restart.'+str(n)+'.nc') 
  39 
  40 
  41 mask=[a.GetTag()==1 for a in atomslist[0]]
  42 band = NudgedElasticBand([Subset(atoms, mask=mask) for atoms in atomslist])
  43 
  44 # Create a quickmin object:
  45 relax = QuasiNewton(band,fmax=0.02,logfilename='restart.log') 
  46 # relax.ReadHessian('hessian.pickle')
  47 
  48 # Create a trajectory for the each image
  49 listoftraj = []
  50 for n in range(len(atomslist)): 
  51     path = NetCDFTrajectory('image-restart'+str(n)+'.nc', atomslist[n])
  52     listoftraj.append(path) 
  53     path.Update()
  54     relax.Attach(path) 
  55 
  56 
  57 relax.Converge()
  58 
  59 # make a trajectory of the final configurations in neb-path.nc
  60 atomslist = [listoftraj[n].GetListOfAtoms(-1) for n in range(len(listoftraj))]
  61 
  62 atom0 = atomslist[0]
  63 newtraj = NetCDFTrajectory('neb-path-restart.nc',atom0) 
  64 newtraj.Update()
  65 
  66 for atom in atomslist[1:]: 
  67     atom0.SetCartesianPositions(atom.GetCartesianPositions())
  68     atom0.SetCalculator(atom.GetCalculator())
  69     newtraj.Update()