Calculation of migration barriers using NEB method

In this example we show how to calculate migration barrier on the example of LiFePO4

Important notes

  • Creation of intermediate images on server is done with nebmake.pl script. By using add_neb(upload_vts = 1) it will be uploaded to ~/tools/vts/

  • initialize header.PATH2NEBMAKE constant with path to nebmake.pl; either in header.py or project_conf.py

  • To use void search you should compile findpores.cpp. For linux use make_findpores script in siman folder. Fow windows use make_findpores.bat (mingw64v5.2 or Visuall C++ should be installed).


In [1]:
import sys
sys.path.append('/home/aksenov/Simulation_wrapper/siman')
from SSHTools import SSHTools
import header
from calc_manage import add_loop, res_loop, smart_structure_read
from database import write_database, read_database
from set_functions import read_vasp_sets
from neb import add_neb
%matplotlib inline

In [2]:
header.ssh_object = SSHTools()
header.ssh_object.setup(user="aksenov",host="10.30.16.62",pkey="/home/aksenov/.ssh/id_rsa")

header.PATH2PROJECT    = 'neb' # path to project relative to your home folder on cluster
header.PATH2POTENTIALS = '/home/aksenov/scientific_projects/PAW_PBE_VASP' #path to VASP POTENTIALS

header.varset['static'].potdir = {3:'Li', 26:'Fe', 15:'P', 8:'O'} #subfolders with required potentials
read_vasp_sets([('ion_relax', 'static', {'ISIF':2, 'IBRION':1, 'NSW':20, 'EDIFFG':-0.025}, 'override')]) #new set 'ion_relax' from 'static' with 'NSW' = 20


Attention! You have chosen to override set ion_relax
 

Warning! You did not change  ISIF  in ion_relax set
 

Warning! You did not change  IBRION  in ion_relax set
 

Out[2]:
{'static': <set_functions.InputSet at 0x7f2620341fd0>,
 'ion_relax': <set_functions.InputSet at 0x7f265c486438>}

In [3]:
st = smart_structure_read(input_geo_file = 'in/mp-19017.POSCAR') #read LiFePO4 structure


-- You provided the following geo file explicitly  in/mp-19017.POSCAR ; Version of file does not matter, I use *curver*= 1 as a new version 

Create Initial and Final configurations and run on server

Function add_neb() creates start and end configurations by removing or adding alkali ion atoms. The start and end configurations have 1 and 2 version numbers, respectively. All intermidiate images have version numbers starting from 3. Before setting run = 1 you can check created configurations in xyz folder: mp-19017_start.xyz and mp-19017_end.xyz

  • it_new and it_new_folder are name and folder for new calculation.
  • ise_new is set with VASP parameters defined above.
  • images is number of images and corenum is number of cores used for calculation; should be dividable by images.
  • upload_vts = 1 upload nebmake.pl to server.

In [4]:
add_neb(st = st, it_new = 'LiFePO4', ise_new = 'ion_relax', it_new_folder = 'LiFePO4/neb',
        images = 3, corenum = 15, upload_vts = 0, run = 0)


-- You have chosen vacancy_creation mode of add_neb tool 

-- Type of atom to move =  Li 

I can suggest you 4 end positions. The distances to them are : [ 0.    3.04  5.74  6.5 ] A
  They are all Li atoms 

-- File xyz/mp-19017_start.xyz was written 

-- File xyz/mp-19017_end.xyz was written 

-- Creating geo files for starting and final configurations (versions 1 and 2)  

-- Warning! File ./LiFePO4/neb/LiFePO4.n3Li1v1/LiFePO4.n3Li1v1.auto_created_starting_position_for_neb_vacancy_creation.1.geo was replaced 

-- Warning! File ./LiFePO4/neb/LiFePO4.n3Li1v1/LiFePO4.n3Li1v1.auto_created_final_position_for_neb_vacancy_creation.2.geo was replaced 

-- Attention, I update number of images in the set to 3 for this calculation;  

-- check_kpoints(): Kpoint   mesh is:  [3, 5, 6] 

-- check_kpoints(): The actual k-spacings are  [ 0.2   0.21  0.22] 


Calculation ('LiFePO4.n3Li1v1', 'ion_relax', 1) successfully created

 

-- check_kpoints(): Kpoint   mesh is:  [3, 5, 6] 

-- check_kpoints(): The actual k-spacings are  [ 0.2   0.21  0.22] 

-- Attention! ngkpt =  [3, 5, 6]  is adopted from struct_des which you provided for it  LiFePO4.n3Li1v1  and kspacing =  0.235 


Calculation ('LiFePO4.n3Li1v1', 'ion_relax', 2) successfully created

 

-- Please be sure that vtsttools are at aksenov@10.30.16.62 /home/aksenov/tools/vts/ 

Out[4]:
'LiFePO4.n3Li1v1'

In [8]:
header.warnings = 'yY' #control warning levels: n, e, y, Y

res_loop('LiFePO4.n3Li1v1', 'ion_relax', [1,2,3,4,5], show = 'fomepp', analys_type = 'neb'  )



Max. F. tot  (meV/A) = 
[1845 2447 1024  738  489  772  231  237  184  126  139  104   55   89   42
   47   37   32   31   27]; 

-- name                          |  energy(eV)|    Vector lenghts (A)   | Stresses (MPa)     | N MD, N SCF    
-- LiFePO4.n3Li1v1.ion_relax.1   | -185.3215  |     10.45;6.09;4.75     |-13095,-13883,-9633 |  20, 7,147    
-- Attention!, SCF was not converged to desirable prec 0.014 > 0.006 meV 



Max. F. tot  (meV/A) = 
[1845 2447 1025  734  492  779  233  240  187  135  151   86   57   82   46
   47   36   38   36   30]; 

-- LiFePO4.n3Li1v1.ion_relax.2   | -185.3211  |     10.45;6.09;4.75     |-13152,-13905,-9671 |  20, 7,149    
-- Directory LiFePO4/neb//LiFePO4.n3Li1v1.ion_relax/01  was created 



Max. F. chain+tot  (meV/A) = 
[1046 1926  794  776  662  593  547  466  379  384  324  348  294  270  248
  226  278  193  182  166]; 

-- LiFePO4.n3Li1v1.ion_relax.3   | -185.1748  |     10.45;6.09;4.75     |-12612,-13980,-8884 |  20, 6,134    
-- Directory LiFePO4/neb//LiFePO4.n3Li1v1.ion_relax/02  was created 



Max. F. chain+tot  (meV/A) = 
[2115 4022 1893 1924 1032  849  773  603  530  850  381  398  340  310  362
  252  249  229  216  198]; 

-- LiFePO4.n3Li1v1.ion_relax.4   | -185.0004  |     10.45;6.09;4.75     |-12148,-14496,-8454 |  20, 7,146    
-- Directory LiFePO4/neb//LiFePO4.n3Li1v1.ion_relax/03  was created 



Max. F. chain+tot  (meV/A) = 
[1052 1932  795  782  663  592  546  465  378  384  326  351  294  271  247
  227  280  190  179  165]; 

-- LiFePO4.n3Li1v1.ion_relax.5   | -185.1748  |     10.45;6.09;4.75     |-12614,-13976,-8877 |  20, 6,138    
[ 0.          1.61416911  3.22825098]
-- plot_mep(): Diffusion barrier = 0.32  eV 

[ 0.          1.61416911  3.22825098]
-- plot_mep(): Diffusion barrier = 0.32  eV 

Out[8]:
({'barrier': 0.32107652044765422}, [])

Start: atom number; End: find interstitial void


In [ ]:
add_neb(RbVl, ise_new = '1u', images = 5, r_impurity = 1.5, search_type = 'existing_voids',
                        i_atom_to_move = 1,  i_void_final = None) #

Start: atom number; End: xred coordinates of void


In [ ]:
# xr_final is coordinates of final position - void
it = add_neb(RbVl, ise_new = '1u', images = 5,
                        i_atom_to_move = 2,  xr_final = [0.2, 0.1, 0.3], i_void_final = 10) # *i_void_final*  is needed here only for producing unique name

Additional help:

  1. produce xyz trajectory calc_manage.py
     res_loop(show = 'neb_geo', params = {'mep_shift_vector'}) 
         neb_analysis()
             write_xyz(sts = sts) # write traectory