View source for Shuffle-Glide dislocation MD and NEB

<P ALIGN="CENTER">
<FONT SIZE="+3" color="darkred"><STRONG>
Shuffle Glide Dislocation Complex: NEB and MD </STRONG></font></P>
<DIV>
<P ALIGN="CENTER"><STRONG>Xiaohan Zhang and Wei Cai</STRONG></P>
</DIV>
<P ALIGN="CENTER"> Created Mar, 2018, Last modified Mar, 2018</P>
<P>

This tutorial corresponds to paper ``Shuffle-Glide Dislocation Complex Nucleation in Silicon Thin Film''. It explains how to calibrate and specify SW and MEAM-lammps potential parameters in both MD++ and LAMMPs. It explains how to set up and run NEB calculations to measure the energy barrier of a shuffle-glide dislocation complex nucleated in a thin film with a surface pit. Finite temperature  MD simulations are performed in Lammps to capture nucleation events as validation to energy barrier calculations. 

<HR>

==Potential for Pure Elements==

===SW Potential for Si===
We calibrate both SW (1985) and SW (1992) potentials in both MD++ and Lammps. The modified version of SW fits to cohesive energy. Their parameters are listed as follows:

MD++:
/* original Si version PRB 31, 5262 (1985) */
       aa=15.27991323; bb=11.60319228; plam=45.51575;
       pgam=2.51412; acut=3.77118; pss=2.0951; rho=4.0;

/* modified Si parameters PRB 46, 2250 (1992) */
       aa=16.31972277; bb=11.60319228; plam=48.61499998;
       pgam=2.51412;  acut=3.77118; pss=2.0951; rho=4.;

Lammps:

/* sw, Stillinger and Weber,  Phys. Rev. B, v. 31, p. 5262, (1985) */

         Si Si Si 2.1674166667   2.0951  1.80  21.0  1.20  -0.333333333
         7.049827  0.602225  4.0  0.0 0.0

/* PRB 46, 2250 (1992) */

         Si Si Si 2.1683  2.0951  1.80  22.4207904718  1.20  -0.333333333333
         7.5265059125  0.6022245574  4.0  0.0 0.0

===MEAM Potential for Si===

We use the 'Siz' potential as those used in Kang, et al "Size and Temperature Effects on Brittle and Ductile Fracture of Silicon Nanowires", International Journal of Plasticity, 26, 1387 (2010" and "Brittle and Ductile Fracture of Semiconductor Nanowires – Molecular Dynamics Simulations", Philosophical Magazine, 87, 2169, (2007)." 


MD++:

The main parameters in the MEAM potential is specified in the '''meamf''' file. (In MD++, this file is in the potentials/MEAMDATA folder.)  The lines correspond to 'Siz' is given below.  

       'Siz'        'dia'   4.      14           28.086
       4.87        4.4     5.5     5.5          5.5   5.431   4.63    1.
       1.0         3.13            4.47         -1.80         1.600   0

corresponding to the data format of:
      elt        lat     z       ielement     atwt
      alpha      b0      b1      b2           b3    alat    esub    asub
      t0         t1              t2           t3            rozero  ibar

Note that  the nearest neighbor distance <math> R_i^0 </math> = '''alat'''  <math>\times \sqrt{3}/4</math> for the diamond cubic structure.

<math>\rho_0^{\rm Si}</math> = '''rozero''' will be important only for cross-potential. And note that this is the only different from Si4 line. 

'''ibar''' is a setting used in the equation of state (EOS), which is explained in SiGe potential tutorial (http://micro.stanford.edu/wiki/MEAM_Potential_for_Si-Ge)



Lammps:

      'Siz'        'dia'   4.      14           28.086
      4.87        4.4     5.5     5.5          5.5   5.431   4.63    1.
      1.0         3.13            4.47         -1.80         1.600   0

corresponds to the data format of 

      DATE: 2012-06-29 DATE: 2007-06-11 CONTRIBUTOR: 
      Greg Wagner, gjwagne@sandia.gov 
      CITATION: Baskes, Phys Rev B, 46, 2727-2742 (1992)
      meam data from vax files fcc,bcc,dia    11/4/92
      elt        lat     z       ielement     atwt
      alpha      b0      b1      b2           b3    alat    esub    asub
      t0         t1              t2           t3            rozero  ibar

===Test case 1: Potential Calibration for SW and MEAM===
This test case is developed in MD++, but also use lammps executable, to evaluate potential energy of a bulk system and a surface pit structure with both SW (1985, 1992) and MEAM potential. To run the test case, follow the steps below,

Take mc2 as an example, one module load the following:

       ::> module load mvapich2/2.0rc1-intel-14, intel/14
       ::> cd MD++.git (svn)
       ::> make sworig build=R SYS=mc2_mpich
       ::> make sw build=R SYS=mc2_mpich
       ::> make meam-lammps build=R SYS=mc2
       ::> bin1/sworig_mc2_mpich scripts/work/sw_meam_calibrations/disl_nuc_hetero.tcl  1 0 1000   

 (0 and 1000 are useless, but need to be there otherwise the script will complain missing parameters)

       ::> bin1/sw_mc2_mpich scripts/work/sw_meam_calibrations/disl_nuc_hetero.tcl  1 0 1000  
       ::> meam-lammps_mc2 scripts/work/sw_meam_calibrations/disl_nuc_hetero.tcl 1 0 1000

change in the script the configuration to read in, choosing from runs/sw_meam_calibrations/bulk_sw.cn, SiGeHole_sw.cn. (There are also bulk_meam.cn, SiGeHole_meam.cn generated from disl_nuc_hetero.tcl with meam-lammps)

From MD++, the numbers should be exactly the following:

       MEAM + bulk_sw.cn :            -1.024081.9
       MEAM + SiGeHole_sw.cn:      -9.97015.8
       SWorig + bulk_sw.cn :           -9.58832.6
       SWorig + SiGeHole_sw.cn :  -9.32569.7
       SW       +  bulk_sw.cn :         -1.024081.9
       SW       + SiGeHole_sw.cn :  -9.96031.739


From lammps,  

       MEAM + bulk_sw.cn :              -1024081.9
       MEAM + SiGeHole_sw.cn:       -997015.9
       SWorig + bulk_sw.cn :              -958832.17 
       SWorig + SiGeHole_sw.cn :      -932569.26
       SW       +  bulk_sw.cn :              -1024081.9
       SW       + SiGeHole_sw.cn :      -996031.74

==NEB calculation work flow==

      ::> cd MD++.git
      ::> module list
      Currently Loaded Modulefiles:
      1) null                       2) intel/14                   3) mvapich2/2.0rc1-intel-14
      ::> make sworig build=R SYS=mc2_mpich

      ::> bin1/sworig_mc2_mpich scripts/work/he-si-tf-pit-neb/he-si-tf-pit-special-3-sworig/disl_nuc_hetero.tcl 0 0 001 1
      ::> bin1/sworig_mc2_mpich scripts/work/he-si-tf-pit-neb/he-si-tf-pit-special-3-sworig/disl_nuc_hetero.tcl 1 0.0520 001 1
      ::> bin1/sworig_mc2_mpich scripts/work/he-si-tf-pit-neb/he-si-tf-pit-special-3-sworig/disl_nuc_hetero.tcl 4 0.0520 001 1

The first line executes a thread that 
1) creates a thin film of (001) free surface
2) reconstruct surface atoms by perturb the atom positions and forming dimers
3) apply along [110] direction a compression strain by the increments of 0.001 from 0 to 0.06, i.e., 6% compressive strain. 

The second line reads in initial config of (3) from above and create a shuffle-glide dislocation complex beneath top free surface, which is then relaxed for several steps. The output configuration can be examined by importing w-loop-compression-surf001-eps0.0520.cfg to ovito/DXA. 

The third line needs to be launched within pbs file `cause it execute sworig_mc2_mpich in MPI mode, calling stringrelax_parallel (). The energy barrier curve can be  visualized by using octave/matlab such as,

===Test case 2: Energy barrier with SW for strain 5.2%===

       ::> load runs/he-si-tf-pit-neb/he-si-tf-pit-special-3-sworig/w-disl-nuc-hetero-001-0.0520.24/stringeng.out
       ::> plot(stringeng(1:100:end, 3:5:end)', '*')

Apply this work flow for different strain, one gets the energy barrier as a function of applied strain.

Note that, dislocation loop size of status == 1 can be adjusted manually in the script by specifying w in make_glide_dislocation_loop_1, for example, minor adjustment is needed to get a smooth energy barrier - strain curve.


==Finite temperature MD with Lammps==

We use lammps-30Jul16 version. 

     ::> cd lammps-30Jul16/lib/meam
     ::> make -f Makefile.gfortran (for meam)
     ::> cd lammps-30Jul16/src
     ::> make yes-manybody
     ::> make yes-meam
     ::> make mpi 


First, one estimates thermal expansion ratios for a bulk system of si using MEAM and sworig. 

     ::> cd MD++.git/
     ::> bin1/sworig_mc2 scripts/work/he-si-tf-pit-md/thermal_init_config/therm_exp_estimates/bulk_generator.tcl 0
     ::> cd /home/xzhang11/Planet/Libs/MD++.git/scripts/work/he-si-tf-pit- 
          md/thermal_init_config/therm_exp_estimates/thermal_equilibrate_load_hold_si_800K
     ::> lmp_mpi < thermal_equilibrate.sw (meam)

From initial and final equilibrated states, one gets thermal expansion. And one records the coefficents in 
         scripts/work/he-si-tf-pit-md/thermal_init_config/therm_initconfig_generator.tcl

     ::> sw_mc2  scripts/work/he-si-tf-pit-md/thermal_init_config/therm_initconfig_generator.tcl 0 1 1000 

To run the lammps scripts, first create surface pit structure with MD++ as follows,

     ::> cd MD++.git/scripts/work/he-si-tf-pit-md/he-si-pit-800K-54-thermal-relax-nvt
     ::> mpirun -np $ncpu lammps-30Jul16/src/lmp_mpi < input.meam

===Test case 3: Nucleation event with SW for strain 5.2%===

       ::> load runs/he-si-tf-pit-neb/he-si-tf-pit-special-3-sworig/w-disl-nuc-hetero-001-0.0520.24/stringeng.out
       ::> plot(stringeng(1:100:end, 3:5:end)', '*')

Apply this work flow for different strain, one gets the energy barrier as a function of applied strain.