Liquid Si - Freezing: Difference between revisions
| Line 145: | Line 145: | ||
=== Diffusion === | === Diffusion === | ||
The diffusion coefficient is given as | |||
<math> D=\frac{\langle x^{2} \rangle} {2 t}</math> | |||
where t defines time and <math>\langle x^{2} \rangle</math>. In our case we calculate the above equation as follows | |||
<math> D=\frac{\langle [x(t)-x(0)]^{2} \rangle}{2 \Delta t} </math>. | |||
To analyse the diffusion behaviour at a certain temperature T, the data read from {{TAG|XDATCAR}}.T can be processed | To analyse the diffusion behaviour at a certain temperature T, the data read from {{TAG|XDATCAR}}.T can be processed | ||
Revision as of 12:39, 3 July 2019
Overview >Liquid Si - Standard MD > Liquid Si - Freezing > Nucleophile Substitution CH3Cl - Standard MD > Nuclephile Substitution CH3Cl - mMD1 > Nuclephile Substitution CH3Cl - mMD2 > Nuclephile Substitution CH3Cl - mMD3 > Nuclephile Substitution CH3Cl - SG > Nuclephile Substitution CH3Cl - BM > List of tutorials
Task
In this example the goal is to simulate the freezing of liquid Si.
Input
POSCAR
Si
15.12409564534287297131
0.5000000000000000 0.5000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.5000000000000000
48
Direct
0.8550657259653851 0.3204575801875221 0.6180363868822553
0.6045454476433229 0.0546379652195404 0.1629680405553871
0.4803889256776521 0.2999635319377835 0.0131251454718051
0.8413504226620471 0.7598095803296524 0.1917781560970181
0.9754163118144437 0.6134171268457649 0.7421364242876367
0.2668229391055025 0.0066502741664650 0.0031140604380929
0.8935777664000575 0.3324172908647429 0.9535738516718881
0.0527608886321274 0.5249316429131962 0.5293744880144071
0.4396089233132741 0.7564833235979471 0.5665855438788387
0.5907859878830199 0.5198033580597228 0.3581725847640679
0.2120832721474721 0.4042899613004446 0.7921535013319151
0.0225803885096466 0.8414911198321031 0.1209255489569852
0.0992500701525566 0.3917384466892963 0.3612433325214984
0.9673794138223195 0.5206425706394114 0.1719623236201897
0.2774602656926126 0.8480860088162007 0.2673309412777037
0.0196991774214161 0.8282178425383616 0.6986213756952502
0.3570927152895376 0.2951488295546784 0.2651851032568589
0.1663829731894614 0.9766237917413699 0.6051764245375237
0.4931841331696695 0.8689890620771937 0.2612357008392290
0.8006473407426477 0.1033419073227807 0.4706563716777467
0.0161340851939779 0.9953827418297991 0.8853439845676159
0.7827740166661069 0.1821830067208054 0.9399555168314748
0.0720651739141343 0.2539424963694544 0.6857919074323433
0.4443385370769313 0.0486404637002326 0.4180706114402839
0.7055263679666055 0.6802623819082319 0.7983614866719116
0.2237125282521105 0.4055474352416297 0.0077044950891134
0.2963682069847125 0.5771265542042112 0.2019757061665083
0.2782449529809642 0.0451513130915826 0.7644934848784113
0.9312079203181675 0.9090938018377080 0.3429249881187518
0.6341882597200124 0.2969253226419481 0.3227590981305088
0.3587691103780569 0.1061057273904179 0.0931868777500710
0.8710437838676732 0.6541301230631744 0.4261617089364881
0.6784300588817769 0.3263889355408940 0.5560491395978739
0.5597052314845080 0.0174390112509929 0.6129003207931863
0.0595962318875451 0.1019295953521402 0.3340999072062676
0.7689671766774326 0.1768870209149794 0.1604177484299765
0.9603661624482890 0.3311649224573259 0.1439224909303592
0.3792868784787023 0.2806150985211180 0.4921541531665999
0.8079860889823454 0.9194188799048340 0.9131036494263627
0.3002081239026374 0.7834053620019006 0.8650323716139056
0.4704528574512951 0.7221628305989689 0.9746107190983403
0.2886552568292480 0.5927625600330780 0.4239421203107919
0.4116743942942291 0.2198943758058664 0.7072597030225044
0.2104494234814825 0.6457654201409418 0.8275863924787099
0.6784628197745537 0.7205455185203838 0.1093053357228383
0.6344130299021448 0.1650970001101275 0.8037018707797643
0.3965793440603315 0.5364088146415013 0.6064549771969059
0.6686412136025504 0.7848666926903073 0.5681234351534038
INCAR
SYSTEM = Si # electronic degrees LREAL = A # real space projection PREC = Normal # chose Low only after tests EDIFF = 1E-5 # do not use default (too large drift) ISMEAR = -1 ; SIGMA = 0.130 # Fermi smearing: 1500 K 0.086 10-3 ALGO = Very Fast # recommended for MD (fall back ALGO = Fast) MAXMIX = 40 # reuse mixer from one MD step to next ISYM = 0 # no symmetry NELMIN = 4 # minimum 4 steps per time step, avoid breaking after 2 steps # MD (do little writing to save disc space) IBRION = 0 # main molecular dynamics tag NSW = 400 # number of MD steps POTIM = 3 # time step of MD NWRITE = 0 # controls output NBLOCK = 10 # after ten steps pair correlation function is written out LCHARG = .FALSE. # no charge density written out LWAVE = .FALSE. # no wave function coefficients written out TEBEG = $i # starting temperature for MD TEEND = $i # end temperature for MD # canonic (Nose) MD with XDATCAR updated every 10 steps MDALGO = 2 ä switch to select thermostat SMASS = 3 # Nose mass ISIF = 2 # this tag selects the ensemble in combination with the thermostat
- Most of the tags here are very similar to the tags used in the previous example (Liquid Si - Standard MD.
- A stepwise cooling will be applied in this example via a script where $i for TEBEG and TEEND will be replaced in each calculation (see below).
KPOINTS
test 0 0 0 monk 1 1 1 0 0 0
- A single k-point is sufficient in this example.
Calculation
We will execute the cooling stepwise so several calculations at different temperatures are required in this calculation. The INCAR is created with a script for each temperature and run separately. After each step the important files are saved to file.$i, where $i are the temperatures ranging from 2000 to 800 K in steps of 100 K. The script running the calculations looks like the following:
for i in 2000 1900 1800 1700 1600 1500 1400 1300 1200 1100 1000 900 800 do cat >INCAR <<! SYSTEM = Si # electronic degrees LREAL = A # real space projection PREC = Normal # chose Low only after tests EDIFF = 1E-5 # do not use default (too large drift) ISMEAR = -1 ; SIGMA = 0.130 # Fermi smearing: 1500 K 0.086 10-3 ALGO = Very Fast # recommended for MD (fall back ALGO = Fast) MAXMIX = 40 # reuse mixer from one MD step to next ISYM = 0 # no symmetry NELMIN = 4 # minimum 4 steps per time step, avoid breaking after 2 steps # MD (do little writing to save disc space) IBRION = 0 # main molecular dynamics tag NSW = 400 # number of MD steps POTIM = 3 # time step of MD NWRITE = 0 # controls output NBLOCK = 10 # after ten steps pair correlation function is written out LCHARG = .FALSE. # no charge density written out LWAVE = .FALSE. # no wave function coefficients written out TEBEG = $i # starting temperature for MD TEEND = $i # end temperature for MD # canonic (Nose) MD with XDATCAR updated every 10 steps MDALGO = 2 ä switch to select thermostat SMASS = 3 # Nose mass ISIF = 2 # this tag selects the ensemble in combination with the thermostat ! mpirun -np 2 /path/to/your/vasp/executable cp XDATCAR XDATCAR.$i cp OUTCAR OUTCAR.$i cp PCDAT PCDAT.$i cp CONTCAR CONTCAR.$i cp POSCAR POSCAR.$i cp OSZICAR OSZICAR.$i cp CONTCAR POSCAR done
- Before running the script one has to replace "'/path/to/your/vasp/executable'" by the path to his "'vasp_gam'" executable. The script is then simply starte by typing "'./script'" on the command line.
Diffusion
The diffusion coefficient is given as
where t defines time and . In our case we calculate the above equation as follows
.
To analyse the diffusion behaviour at a certain temperature T, the data read from XDATCAR.T can be processed using the script diffusion (caution - this script is written for VASP version 4.x output):
awk <XDATCAR >diffusion.xy '
#
# simple module function
#
function mod(x,y) { return x-int(x/y)*y }
function minim(x) { return mod(x+2.5,1.0)-0.5 }
#
# calculate mean square displacement
#
function diff() {
d=0
for (ion=1; ion<=ions; ion++) {
dx=minim(xn[ion]-x[ion])
dy=minim(yn[ion]-y[ion])
dz=minim(zn[ion]-z[ion])
xn[ion]=x[ion]+dx
yn[ion]=y[ion]+dy
zn[ion]=z[ion]+dz
d=d+(xn[ion]-x0[ion])*(xn[ion]-x0[ion])*a1*a1
d=d+(yn[ion]-y0[ion])*(yn[ion]-y0[ion])*a2*a2
d=d+(zn[ion]-z0[ion])*(zn[ion]-z0[ion])*a3*a3
}
# d=d/(set*t)/6
d=d/6
print set*t,d
}
#
# set the number of ions
#
NR==1 { ions = $1 }
NR==2 { a1=$2*10^10 ; a2=$3*10^10 ; a3=$4*10^10 ; t=$5*10^12 }
#
# at this point a complete set of ionic positions has been found
#
mod(NR-6,ions+1)==0 {
if (set>=2) diff()
if (set==1) {
for (ion=1; ion<=ions; ion++) {
x0[ion]=xn[ion]
y0[ion]=yn[ion]
z0[ion]=zn[ion]
}
}
for (ion=1; ion<=ions; ion++) {
x[ion]=xn[ion]
y[ion]=yn[ion]
z[ion]=zn[ion]
}
head=headn
headn=$0
set=set+1
}
# store coordinates
mod(NR-6,ions+1)>0 {
ion=mod(NR-6,ions+1)
xn[ion]=$1
yn[ion]=$2
zn[ion]=$3
}
'
Pair correlation function
The pair-correlation function provides information about the probability of finding two atoms at a given distance . The pair-correlation function written on PCDAT.T should be processed using the script PCDATtoPCDATxy:
awk <PCDAT >PCDAT.xy '
NR==8 { pcskal=$1}
NR==9 { pcfein=$1}
NR>=13 {
line=line+1
if (line==257) {
print " "
line=0
}
else
print (line-0.5)*pcfein/pcskal,$1
}
'
Mind: You will have to set the correct path to your VASP executable and invoke VASP with the correct command (e.g., in the above: mpirun -np 2).