Liquid Si - Standard MD: Difference between revisions
(→150 ps) |
(→150 ps) |
||
Line 156: | Line 156: | ||
gnuplot -e "set terminal jpeg; set xlabel 'r(Ang)'; set ylabel 'PCF'; set style data lines; plot 'PCDAT.150ps'" > PC_150ps.jpg | gnuplot -e "set terminal jpeg; set xlabel 'r(Ang)'; set ylabel 'PCF'; set style data lines; plot 'PCDAT.150ps'" > PC_150ps.jpg | ||
The obtained pair correlation function should look like in Fig. 1. Solids usually show sharp peaks in the pair correlation function since the ions only vibrate around fixed positions in the crystal lattice. In the liquid or amorphous state the distances are much more diffuse and one usually would expect no far order (but both can have near order). We see that many sharp peaks arise in the pair correlation function at higher distances, so that after 150 ps we have not molten the material. | The obtained pair correlation function should look like in Fig. 1. Solids usually show sharp peaks in the pair correlation function since the ions only vibrate around fixed positions in the crystal lattice. In the liquid or amorphous state the distances are much more diffuse and one usually would expect no far order (but both can have near order). We see that many sharp peaks arise in the pair correlation function at higher distances, so that after 150 ps we have not molten the material. | ||
[[File:Energy_150ps.jpg|300px|thumb|Fig. 2: Total energy vs number of steps in first 150 ps.]] | [[File:Energy_150ps.jpg|300px|thumb|Fig. 2: Total energy vs number of steps in first 150 ps.]] |
Revision as of 14:34, 21 June 2019
Task
Generating liquid Si by melting of the crystalline structure via molecular dynamics.
Input
POSCAR
- We start this example by making a POSCAR using the conventional unit cell with 8 atoms which should look like this:
Si cubic diamond conventional cell 5.43100000000000 1.00000000 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 1.00000000 Si 8 Direct 0.00000000 0.00000000 0.00000000 0.75000000 0.25000000 0.75000000 0.00000000 0.50000000 0.50000000 0.75000000 0.75000000 0.25000000 0.50000000 0.00000000 0.50000000 0.25000000 0.25000000 0.25000000 0.50000000 0.50000000 0.00000000 0.25000000 0.75000000 0.75000000
- We obtain a sufficiently large supercell (2x2x2 containing 64 atoms) by following the description in: Preparing a Super Cell.
- The new POSCAR file of the two 2x2x2 super cell of the conventional cell should look like this:
Si cubic diamond 2x2x2 super cell of conventional cell 5.43090000000000 2.00000000 0.00000000 0.00000000 0.00000000 2.00000000 0.00000000 0.00000000 0.00000000 2.00000000 Si 64 Direct 0.00000000 0.00000000 0.00000000 0.50000000 0.00000000 0.00000000 0.00000000 0.50000000 0.00000000 0.50000000 0.50000000 0.00000000 0.00000000 0.00000000 0.50000000 0.50000000 0.00000000 0.50000000 0.00000000 0.50000000 0.50000000 0.50000000 0.50000000 0.50000000 0.37500000 0.12500000 0.37500000 0.87500000 0.12500000 0.37500000 0.37500000 0.62500000 0.37500000 0.87500000 0.62500000 0.37500000 0.37500000 0.12500000 0.87500000 0.87500000 0.12500000 0.87500000 0.37500000 0.62500000 0.87500000 0.87500000 0.62500000 0.87500000 0.00000000 0.25000000 0.25000000 0.50000000 0.25000000 0.25000000 0.00000000 0.75000000 0.25000000 0.50000000 0.75000000 0.25000000 0.00000000 0.25000000 0.75000000 0.50000000 0.25000000 0.75000000 0.00000000 0.75000000 0.75000000 0.50000000 0.75000000 0.75000000 0.37500000 0.37500000 0.12500000 0.87500000 0.37500000 0.12500000 0.37500000 0.87500000 0.12500000 0.87500000 0.87500000 0.12500000 0.37500000 0.37500000 0.62500000 0.87500000 0.37500000 0.62500000 0.37500000 0.87500000 0.62500000 0.87500000 0.87500000 0.62500000 0.25000000 0.00000000 0.25000000 0.75000000 0.00000000 0.25000000 0.25000000 0.50000000 0.25000000 0.75000000 0.50000000 0.25000000 0.25000000 0.00000000 0.75000000 0.75000000 0.00000000 0.75000000 0.25000000 0.50000000 0.75000000 0.75000000 0.50000000 0.75000000 0.12500000 0.12500000 0.12500000 0.62500000 0.12500000 0.12500000 0.12500000 0.62500000 0.12500000 0.62500000 0.62500000 0.12500000 0.12500000 0.12500000 0.62500000 0.62500000 0.12500000 0.62500000 0.12500000 0.62500000 0.62500000 0.62500000 0.62500000 0.62500000 0.25000000 0.25000000 0.00000000 0.75000000 0.25000000 0.00000000 0.25000000 0.75000000 0.00000000 0.75000000 0.75000000 0.00000000 0.25000000 0.25000000 0.50000000 0.75000000 0.25000000 0.50000000 0.25000000 0.75000000 0.50000000 0.75000000 0.75000000 0.50000000 0.12500000 0.37500000 0.37500000 0.62500000 0.37500000 0.37500000 0.12500000 0.87500000 0.37500000 0.62500000 0.87500000 0.37500000 0.12500000 0.37500000 0.87500000 0.62500000 0.37500000 0.87500000 0.12500000 0.87500000 0.87500000 0.62500000 0.87500000 0.87500000
KPOINTS
K-Points 0 Gamma 1 1 1 0 0 0
- Since a sufficiently large super cell is used in this example, it is ok in this case to use only a single k-point in the calculations. Hence it is also possible to use the -point only version which is significantly faster than the standard version.
INCAR
ISMEAR = 0 IBRION = 0 MDALGO = 2 ISIF = 2 SMASS = 1.0 SIGMA = 0.1 LREAL = Auto ALGO = VeryFast PREC = Low ISYM = 0 TEBEG = 2000 NSW = 50 POTIM = 3.0 NCORE = 2
- To select a molecular dynamics calculation set IBRION=0.
- By selecting MDALGO=2 and ISIF=2 we select the NVT ensemble using the Nose-Hoover thermostat.
- The tag SMASS specifies the Nose mass, which is a ficitional mass for the fictional coordinate of the heat bath. The choice of SMASS=1.0 should work well for this tutorial.
- Since we are dealing with a super cell, we set LREAL=Auto. In this mode the projection operators are evaluated in real space. This should speed up the calculation while being slightly less accurate then the evaluation of the operators in reciprocal space.
- To significantly speed up the calculations we use ALGO=VeryFast and PREC=Low. This is perfectly ok for this tutorial example but for more precise results these flags should be used with caution.
- A time step of 3 femtoseconds (POTIM=3.0) is employed in this example, which should be ok for many applications of Si.
- The tag NCORE=2 specifies that the parallelization is done such that 2 cores share the work on one orbital. This means that for e.g. 8 cores 4 different orbitals would be treated simultaneously, where for each orbital two plane-wave coefficients would be calculated simultaneously.
Calculation
The calculation is started from the perfect crystal. Since the chosen temperature at 2000K is significantly above the known melting temperature at around 1400 K the melting should be achieved relatively quickly.
It is suggested to run this calculation using the -point only version, since we have only one k point. Then it should be a very quick calculation on eight nodes:
mpirun -np 8 vasp_path/vasp_gam
When the solid melts the crystal structure and the distances between the atoms change. This can be well monitored by looking at the pair correlation function (or radial distribution function). We will also monitor the energy conservation to see if we are well equilibrated.
150 ps
First we run the calculation for 150 ps.
![](/wiki/images/thumb/e/e5/PC_150fs.jpeg/300px-PC_150fs.jpeg)
- Pair correlation function:
The pair correlation function is written out to the PCDAT file. The abscissa of that file is within mesh points of a selected grid and need to be converted to . This is done by invoking the following short awk script on the command line:
awk <PCDAT >PCDAT.150ps ' NR==8 {pcskal=$1} NR==9 {pcfein=$1} NR>=13 {line=line+1; print (line-0.5)*pcfein/pcskal,$1} '
This produces the file PCDAT.150ps which contains the pair correlation function against the radius in Angstrom after the first 150 ps (50 timesteps NSW=50 with a stepsize of 3 ps POTIM=3). Now we will plot the pair correlation function in the gnuplot window:
gnuplot -e "set xlabel 'r(Ang)'; set ylabel 'PCF'; set style data lines; plot 'PCDAT.150ps'; pause -1"
or to a file (in the remainder we just show how to plot to a file):
gnuplot -e "set terminal jpeg; set xlabel 'r(Ang)'; set ylabel 'PCF'; set style data lines; plot 'PCDAT.150ps'" > PC_150ps.jpg
The obtained pair correlation function should look like in Fig. 1. Solids usually show sharp peaks in the pair correlation function since the ions only vibrate around fixed positions in the crystal lattice. In the liquid or amorphous state the distances are much more diffuse and one usually would expect no far order (but both can have near order). We see that many sharp peaks arise in the pair correlation function at higher distances, so that after 150 ps we have not molten the material.
- Energy conservation:
We will also output the total energy for each molecular dynamics step by invoking the command:
grep "free energy" OUTCAR|awk ' {print $5}' > energy.dat
We will now plot the energy using gnuplot (the user can of course use the plotting program of choice):
gnuplot -e "set terminal jpeg; set style data lines; plot 'energy.dat'" > energy.jpg
300 ps
We repeat the calculation for another 150 ps.
Before we run the calculation we need to copy the new positions and velocities in CONTCAR to POSCAR. Then rerun the calculation using the same command as above. The pair correlation function is written to the file PCDAT.300ps:
awk <PCDAT >PCDAT.300ps ' NR==8 {pcskal=$1} NR==9 {pcfein=$1} NR>=13 {line=line+1; print (line-0.5)*pcfein/pcskal,$1} '
The new energies are concatenated to the "energy.dat" file (mind the ">>" instead of ">" above):
grep "free energy" OUTCAR|awk ' {print $5}' > energy.dat