Reconstruction of surface with MC

In this tutorial we show how to make DFT Monte-Carlo for surfaces. The Monte Carlo method creates anti-site defects near surface.

In order to use MC:

  1. Install Siman version >0.9.5

sudo python -m pip install siman

  1. Copy installed siman *.py soruces to cluster ~/tools/siman

In [1]:
from siman import header
from siman.calc_manage   import add, res, smart_structure_read
from siman.geo import create_surface2, ortho_vec, create_supercell
from siman.header import db
from siman.analysis import suf_en
from siman.classes import CalculationVasp
from siman.set_functions import read_vasp_sets
header.PATH2POTENTIALS = '/home/aksenov/scientific_projects/PAW_PBE_VASP'

Reading  /home/aksenov/

1. Read bulk structure and create new set

In [2]:
st = smart_structure_read('NaNiO2/POSCAR_NaNiO2_c2_m_su_4uiAg_100_end') # read bulk NaNiO2 with c2/m structure
sconv = st.get_conventional_cell()

dftu_packet = {'ISTART'   :1,   'ICHARG':1,  'LDAUTYPE':2, 'LASPH':'.TRUE.', 
                'LDAUPRINT':2, 'LMAXMIX' :4, 'LDAU' :'.TRUE.',
                'LDAUL':{'Ti':2,   'Co':2  , 'Fe':2  , 'Ni':2  , 'Mn':2  , 'V':2   , 'Cr':2 },
                'LDAUU':{'Ti':0,   'Co':3.4, 'Fe':4.0, 'Ni':6.2, 'Mn':3.9, 'V':3.1 , 'Cr':3.5, 'Fe/S':1.9 },
                'LDAUJ':{'Ti':0.0, 'Co':0.0, 'Fe':0.0, 'Ni':0.0, 'Mn':0.0, 'V':0.0 , 'Cr':0.0, 'Fe/S':0   } }

mag_packet = {
    'GGA_COMPAT': '.FALSE.',
    'LORBIT':11, #more info
    'magnetic_moments':{'Ti':0.6, 'V':5, 'Fe':5, 'Co':5, 'Mn':5, 'Ni':5, 'Cr':5 }

calc_pack = {'KSPACING':0.3, 'ENCUT':400, 'ENAUG':400*1.75,  }

surface_pack4 = {'ISYM':0, 'PREC':'Accurate', 'NSW':50, 'NELM':100, 'EDIFF':1e-04, 'ISMEAR':0, 'SIGMA':0.1, 'IBRION':2, 'EDIFFG':-0.05, 'POTIM':0.2,
                 'AMIN':None, 'AMIX':None, 'BMIX':None, 'NELMIN':4, 'IDIPOL':3, 'LDIPOL':'.TRUE.', 'LVTOT':'.TRUE.', 'ICHARG':1} # 

pack = {}
# print(pack)
read_vasp_sets([('1u6', 'static', pack)])

Warning! You did not change  AMIX  in 1u6 set

Warning! You did not change  AMIN  in 1u6 set

Warning! You did not change  BMIX  in 1u6 set

Warning! You did not change  NELM  in 1u6 set

Warning! You did not change  NELMIN  in 1u6 set

{'1u6': <siman.set_functions.InputSet at 0x7fdf33723cf8>,
 None: <siman.set_functions.InputSet at 0x7fdf33723d30>,
 'opt': <siman.set_functions.InputSet at 0x7fdf43694668>,
 'static': <siman.set_functions.InputSet at 0x7fdf37f0bc18>}

2. Create slab with surface surface

Here the second argument provide the Miller indexes of required surface

 - 'shift' is additional shift of input cell along R3
 - 'cut_thickness' allows to remove layers from the top of slab after its construction
 - 'surface_i' is choosing the required termination if more than one exists
 - 'symmetrize' is symmetrizing the slab to find the surface primitive cell 

In [3]:
st_sufNa = create_surface2(sconv, [1, 0, -1], shift = 0.0, 
                           return_one = 1, write_poscar = 0, cut_thickness = 0, 
                           min_vacuum_size = 10, symmetrize = 0, 
                           min_slab_size = 16, surface_i = 0, oxidation = {'Na':'Na+', 'Ni':'Ni3+', 'O':'O2-' })

-- 1 surfaces were generated, choose required surface using *surface_i* argument 

Final structure contains  32 atoms

3. Fix bottom layers of the slab

Here we selective dynamics to fix bottom layers of the slab. This allows to speed up relaxation of atomic positions.
The higlight argument allow to write xyz file to check which atoms would be fixed.

In [4]:
st_sufNaf = st_sufNa.fix_layers(xcart_range = [0, 6.8], highlight = 0)

4. Submit Monte-Carlo job

Choose MC regime with calc_method argument.
The parameters of MC simulation are given in params dic with 'monte' section

- 'normal'   - choose vector normal to surface
- 'thickness' - thickness of slice where Monte-Carlo changes are allowed (from the top surface)
- 'mcsteps' - number of Monte-Carlo steps
- 'temp'    - temperature (K) for Metropolis Algorithm

In [8]:
header.warnings = 'yY'
add('', '1u6', 1, up = 'up2', input_st = st_sufNaf, it_folder = 'NaNiO2/mc1/10-1/', 
    calc_method = 'monte', run = 0,
    params = {'monte':{'normal':2, 'thickness':6, 'mcsteps':20, 'temp':1073}}, show ='fo')

-- Attention!, cluster None is not found, using default CEE 

-- check_kpoints(): Kpoint   mesh is:  [8, 5, 1] 

-- check_kpoints(): The actual k-spacings are  [0.27 0.24 0.2 ] 

-- actualize_set(): Magnetic moments are determined from self.init.magmom: [0.6, 0.6, 5, 5, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 5, 5, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 5, 5, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 5, 5, 0.6, 0.6, 0.6, 0.6] 

-- POSCAR was written to /home/aksenov/Simulation_wrapper/siman/tutorials/NaNiO2/mc1/10-1/// 

-- Attention! ngkpt =  [8, 5, 1]  is adopted from struct_des which you provided for it  and kspacing =  0.3 

Calculation ('', '1u6', 1) successfully created



Understanding results

Go to calculation folder on cluster. You can see two new files monte.log, verbose_log, and ENERGIES
monte.log gives concise information about calculation, while verbose_log contains more details.

ENERGIES file contains number of step and final energy for this step. You can open this file and see if the energy was reduced due to reconstruction

You can see that for each MC step siman creates OUTCAR-n, CONTCAR-n, OSZICAR-n, n-yes/no.pickle. n-yes/no.pickle is serialized CalculationVasp object. If the filename contains yes in its name, than it means that this configuration was accepted with Monte-Carlo algorithm.

To read results use res() as always. To download CONTCAR-n and n-yes/no.pickle files of intermediate steps use 'pickle' value for show argument

In [ ]:
res('', '1u6', 1, show ='pickle')