Revision as of 02:13, 24 November 2008 (view source) Shryu (Talk) (→‎Submit jobs) ← Older edit

Revision as of 20:25, 7 December 2008 (view source) Caiwei (Talk) (uploaded pbs scripts and picture) Newer edit →

---

Revision as of 20:25, 7 December 2008

Computing Melting Point by Free Eneregy Method

Seunghwa Ryu and Wei Cai

This tutorial describes how to compute melting points from free energy methods. We provide the MD++ script files and describe how to use them in detail. The theoretical background is published in

Comparison of Thermal Properties Predicted by Interatomic Potential Models, Modelling and Simulation in Materials Science and Engineering, 16, 085005 (2008). (PDF).

Contents 1 Download files 2 Compile excutable file 3 Submit jobs 4 Analyse data

Download files

First, you need to install MD++ on your computer by following the instructions on MD++ Manuals.

Second, copy files melting_cubic.tcl and melting_noncubic.tcl to your input file directory of MD++. You can do so by

 mkdir -p ~/Codes/MD++/scripts/work/melting
 cd ~/Codes/MD++/scripts/work/melting
 wget http://micro.stanford.edu/mediawiki-1.11.0/images/Melting_cubic.tcl.txt -O melting_cubic.tcl
 wget http://micro.stanford.edu/mediawiki-1.11.0/images/Melting_noncubic.tcl.txt -O melting_noncubic.tcl

Both files contain a big comment section in the beginning that describes in detail the (16) steps to compute the melting point. These steps are fully automated. In this tutorial, we describe how to run these scripts on a parallel cluster (using wcr.stanford.edu as an example) and what kind of results you should expect.

Compile excutable file

For each different potential, MD++ has corresponding excutable files.

Stillinger-Weber Silicon   - sw
Stillinger-Weber Germanium - swge
MEAM                       - meam-lammps, meam-baskes, meam
EAM                        - eam
Finnis-Sinclair            - fs
Tersoff Silicon            - tersoff 

For example, you can compile sworig by doing

 cd ~/Codes/MD++/
 make sworig build=R SYS=wcr

then you will have 'sworig_wcr' file in the directory '~/Codes/MD++/bin'

For other potentials, refer to the explanations in melting_cubic.tcl or email to shryu@stanford.edu for any question.

Submit jobs

Here are the PBS scripts ( melt1.pbs and melt2.pbs)for submitting the jobs on wcr.stanford.edu.

 cd ~/Codes/MD++/
 mkdir runs/Single_Elem_Tm
 wget http://micro.stanford.edu/mediawiki-1.11.0/images/Melt1.pbs.txt -O melt1.pbs
 wget http://micro.stanford.edu/mediawiki-1.11.0/images/Melt2.pbs.txt -O melt2.pbs
 qsub melt1.pbs
 qsub melt2.pbs

Then you'll get free energy data in

~/Codes/MD++/runs/Single_Elem_Tm/SW_Si_Solid/
~/Codes/MD++/runs/Single_Elem_Tm/SW_Si_Liquid/

Analyse data

At the end of the calculation, data files and Matlab files will be automatically generated. Transfer these files to a location where you have Matlab installed. Run the following commands in Matlab and you will get the free energy plot (for solid and liquid phases). The melting point is the temperature at which the two free energy curves cross.

The analysis files are generated in ~/Codes/MD++/runs/Single_Elem_Tm/

Hessian_SW_Si.m
Free_Energy_of_SW_Si.m
tar_SW_Si

to zip necessary data

 cd ~/Codes/MD++/runs/Single_Elem_Tm
 ./tar_SW_Si

then you'll have 'inter_SW_Si.tar'

copy this file to your local computer and unzip. Then, execute

Hessian_SW_Si.m
Free_Energy_of_SW_Si.m

you'll see following output files including all free energy data and error analysis

Solid Free Energy Info
lattice const =    5.43095
strain =    0.00518
E_0 = -2218.74337 eV (   -4.33348 eV/atom)
F_w = -253.34582 eV (   -0.49482 eV/atom)
F_ha = E_0 + F_w = -2472.08919 eV (   -4.82830 eV/atom)
W_re_to_ha =    5.99169 eV (    0.01170 eV/atom)
W_ha_to_re =   -7.30068 eV (   -0.01426 eV/atom)
F_re - F_ha = (W_ha_to_re - W_re_to_ha)/2 =   -6.64619 eV (   -0.01298 eV/atom)
F_re at 1600.000000 = -2478.73538 eV (   -4.84128 eV/atom)
W_T1_to_T11 =  417.87095 eV (    0.81615 eV/atom)
W_T11_to_T1 = -417.87508 eV (   -0.81616 eV/atom)
dW1 = (W_T1_to_T11 - W_T11_to_T1)/2 =  417.87302 eV (    0.81616 eV/atom)
F_re at 2000.000000 = -2605.55266 eV (   -5.08897 eV/atom)

Liquid Free Energy Info
strain =   -0.01894
F_ideal = -919.36940 (   -1.79564 eV/atom)
W_re_to_ga = 1888.21102 eV (    3.68791 eV/atom)
W_ga_to_re = -1887.53306 eV (   -3.68659 eV/atom)
F_re - F_ga = (W_ga_to_re - W_re_to_ga)/2 = -1887.87204 eV (   -3.68725 eV/atom)
W_ga_to_id = -256.29295 eV (   -0.50057 eV/atom)
W_id_to_ga =  256.52177 eV (    0.50102 eV/atom)
F_ga - F_id = (W_id_to_ga - W_ga_to_id)/2 =  256.40736 eV (    0.50080 eV/atom)
F_re at 1800.000000 = -2550.83408 eV (   -4.98210 eV/atom)
W_T2_to_T22 = -583.02785 eV (   -1.13873 eV/atom)
W_T22_to_T2 =  583.06542 eV (    1.13880 eV/atom)
dW2 = (W_T2_to_T22 - W_T22_to_T2)/2 = -583.04664 eV (   -1.13876 eV/atom)
F_re at 1384.615385 = -2386.68534 eV (   -4.66149 eV/atom)

Error Analysis

Solid Side
Err_re_to_ha =    0.02769 eV ( 0.00005407 eV/atom)( 0.628 K)
Err_ha_to_re =    0.02815 eV ( 0.00005498 eV/atom)( 0.638 K)
Err_T1_to_T11 =    0.00374 eV ( 0.00000730 eV/atom)( 0.085 K)
Err_T11_to_T1 =    0.00441 eV ( 0.00000861 eV/atom)( 0.100 K)

Liquid Side
Err_re_to_ga =    0.04793 eV ( 0.00009362 eV/atom)( 1.086 K)
Err_ga_to_re =    0.06601 eV ( 0.00012893 eV/atom)( 1.496 K)
Err_ga_to_id =    0.02557 eV ( 0.00004993 eV/atom)( 0.580 K)
Err_id_to_ga =    0.03421 eV ( 0.00006682 eV/atom)( 0.776 K)
Err_T2_to_T22 =    0.01488 eV ( 0.00002906 eV/atom)( 0.337 K)
Err_T22_to_T2 =    0.01183 eV ( 0.00002311 eV/atom)( 0.268 K)
Max dT =  2.315
Tm= 1695.08000 (K), err= 1.044 Fm =   -4.89779 (eV/atom)
L=      165.8833097832 (eV) Ss=0.308452 Sl=0.406313
L=        0.3239908394 (eV/atom) ss=0.000602 sl=0.000794
L=    31206.7976529696 (J/mol) ss=58.027496 sl=76.437718
L=  1111115.7748689589 (J/kg) ss=2066.064819 sl=2721.559412