VASP Computing Elastic Constants of Si: Difference between revisions
| (2 intermediate revisions by the same user not shown) | |||
| Line 168: | Line 168: | ||
<pre> |
<pre> |
||
#!/bin/bash |
#!/bin/bash |
||
#PBS -N vasp.Si. |
#PBS -N vasp.Si.C44 |
||
#PBS -j oe |
#PBS -j oe |
||
#PBS -l nodes=1:ppn=8,walltime=48:00:00 |
#PBS -l nodes=1:ppn=8,walltime=48:00:00 |
||
| Line 309: | Line 309: | ||
The first line fits the energy data to a quadratic curve and computes the equilibrium lattice constant, cohesive energy and bulk modulus(['''a0, Ecoh, B]''') respectively. The parameter '''2''' specifies the number of atoms in the computation cell and '''8''' is the number of atoms in a conventional unit cell. |
The first line fits the energy data to a quadratic curve and computes the equilibrium lattice constant, cohesive energy and bulk modulus(['''a0, Ecoh, B]''') respectively. The parameter '''2''' specifies the number of atoms in the computation cell and '''8''' is the number of atoms in a conventional unit cell. |
||
The second lines computes the C11 and C12 elastic constants and the third line calculates C44. |
The second lines computes the C11 and C12 elastic constants and the third line calculates C44. Suggested results are C11 = 161 GPa, C12 = 63 GPa and C44 = 76 GPa for calculations with US-LDA PP. |
||
<span style="background:yellow">'''C44 data we got here is much smaller than the paper's result of 105 GPa but is the same as the experimental value of 76 GPa. We suspect maybe 105 GPa is the unrelaxed result since only specifying ISIF=2 doesn't enable relaxation. ''' </span> |
|||
Latest revision as of 05:40, 6 September 2015
VASP: Computing Elastic Constants of Si
Adriano Sanchez，Yanming Wang and Wei Cai
Input files
Here we give an example of how to use VASP to compute the elastic constants C11, C12 and C44 of Si. We performed this calculation on MC2 in PARALLEL mode in ~/Codes/VASP/runs/Si_bulk_C11 and ~/Codes/VASP/runs/Si_bulk_C44 two directories. Each directory contains the following files.
INCAR
PREC = High ISTART = 0 ICHARG = 2 ISMEAR = 1 SIGMA = 0.1 EDIFF = 1E-09 NELM = 40 ENMAX = 500 ENCUT = 500 ISIF = 2 NSW = 100 IBRION = 2
KPOINTS
21x21x21 0 0 = automatic generation of k-points Monkhorst 21 21 21 0 0 0
POSCAR
Here's an example for Si C11 and C12 calculation, which will be changed later in PBS script every iteration.
POSCAR for DC Si 5.402856 accurate equilimrium lattice constant 1.003 0.0 0.0 first Bravais lattice vector 0.0 1.0 0.0 second Bravais lattice vector 0.0 0.0 1.0 third Bravais lattice vector 8 number of atoms per species selective dynamics direct direct or cart (only first letter is significant) 0.00000 0.00000 0.00000 T T T 0.50000 0.50000 0.00000 T T T 0.00000 0.50000 0.50000 T T T 0.50000 0.00000 0.50000 T T T 0.25000 0.25000 0.25000 T T T 0.75000 0.75000 0.25000 T T T 0.25000 0.75000 0.75000 T T T 0.75000 0.25000 0.75000 T T T
Here's an example for Si C44 calculation, which will be changed later in PBS script every iteration.
POSCAR for DC Si 5.402856 accurat equilibrium lattice constant 1.000 1.000 0.0 first Bravais lattice vector -1.0 1.0 0.0 second Bravais lattice vector 0.0 0.0 1.0 third Bravais lattice vector 16 number of atoms per species selective dynamics direct direct or cart (only first letter is significant) 0.0000 0.5000 0.0000 T T T 0.2500 0.7500 0.5000 T T T 0.0000 0.0000 0.0000 T T T 0.5000 0.0000 0.0000 T T T 0.2500 0.2500 0.5000 T T T 0.5000 0.5000 0.0000 T T T 0.7500 0.7500 0.5000 T T T 0.7500 0.2500 0.5000 T T T 0.2500 0.5000 0.2500 T T T 0.5000 0.7500 0.7500 T T T 0.2500 0.0000 0.2500 T T T 0.7500 0.0000 0.2500 T T T 0.5000 0.2500 0.7500 T T T 0.7500 0.5000 0.2500 T T T 1.0000 0.7500 0.7500 T T T 1.0000 0.2500 0.7500 T T T
POTCAR
In this example, we choose the ultrasoft-LDA PP of Si provided by VASP.
Run VASP
To compute the elastic constants, we prepared the following PBS scripts to submit the jobs to MC2 cluster. vasp.pbs.C11 is used to to run vasp repeated with different values for the first component of the first Bravais lattice vector, $a. Here we put 5.402856 as the length of the lattice constant in the script. In generally, you are expected to put a reliable value from your lattice constant calculation with the same PP. The script automatically creates the POSCAR file for every a specified in the script in the range from 0.997 to 1.003. For C44 calculation, the PBS script vasp.pbs.C44 changes the value of $a for the first and second component of the first Bravais lattice vector.
PBS script for C11
#!/bin/bash #PBS -N vasp.Si.C11 #PBS -j oe #PBS -l nodes=1:ppn=8,walltime=48:00:00 #PBS -V ### --------------------------------------- ### BEGINNING OF EXECUTION ### --------------------------------------- echo The master node of this job is `hostname` echo The working directory is `echo $PBS_O_WORKDIR` echo This job runs on the following nodes: echo `cat $PBS_NODEFILE` ncpu=`cat $PBS_NODEFILE | wc -w` echo "Number of processors = $ncpu " ### end of information preamble cd $PBS_O_WORKDIR echo $PWD cd $PBS_O_WORKDIR module load vasp rm WAVECAR for a in 0.997 0.998 0.999 1.000 1.001 1.002 1.003 do cat > POSCAR << FIN POSCAR for DC Si 5.402856 accurate equilibrium lattice constant $a 0.0 0.0 first Bravais lattice vector 0.0 1.0 0.0 second Bravais lattice vector 0.0 0.0 1.0 third Bravais lattice vector 8 number of atoms per species selective dynamics direct direct or cart (only first letter is significant) 0.0000 0.0000 0.0000 T T T 0.5000 0.5000 0.0000 T T T 0.0000 0.5000 0.5000 T T T 0.5000 0.0000 0.5000 T T T 0.2500 0.2500 0.2500 T T T 0.7500 0.7500 0.2500 T T T 0.2500 0.7500 0.7500 T T T 0.7500 0.2500 0.7500 T T T FIN echo "a=$a ncpu=$ncpu" cmd="mpiexec -np $ncpu vasp" $cmd >& vasp.log E=`tail -1 OSZICAR` echo $a $E | sed -s 's/F=//; s/E0=//; s/d E =//;' >> Elatt.C11.dat done
PBS script for C44
#!/bin/bash #PBS -N vasp.Si.C44 #PBS -j oe #PBS -l nodes=1:ppn=8,walltime=48:00:00 #PBS -V ### --------------------------------------- ### BEGINNING OF EXECUTION ### --------------------------------------- echo The master node of this job is `hostname` echo The working directory is `echo $PBS_O_WORKDIR` echo This job runs on the following nodes: echo `cat $PBS_NODEFILE` ncpu=`cat $PBS_NODEFILE | wc -w` echo "Number of processors = $ncpu " ### end of information preamble cd $PBS_O_WORKDIR echo $PWD cd $PBS_O_WORKDIR module load vasp rm WAVECAR for a in 0.997 0.998 0.999 1.000 1.001 1.002 1.003 do cat > POSCAR << FIN POSCAR for DC Si 5.402856 accurate equilibrium lattice constant $a $a 0.0 first lattice vector -1.0 1.0 0.0 second lattice vector 0.0 0.0 1.0 third lattice vector 16 number of atoms per species selective dynamics direct direct or cart (only first letter is significant) 0.0000 0.5000 0.0000 T T T 0.2500 0.7500 0.5000 T T T 0.0000 0.0000 0.0000 T T T 0.5000 0.0000 0.0000 T T T 0.2500 0.2500 0.5000 T T T 0.5000 0.5000 0.0000 T T T 0.7500 0.7500 0.5000 T T T 0.7500 0.2500 0.5000 T T T 0.2500 0.5000 0.2500 T T T 0.0000 0.7500 0.7500 T T T 0.0000 0.2500 0.7500 T T T 0.5000 0.7500 0.7500 T T T 0.2500 0.0000 0.2500 T T T 0.7500 0.0000 0.2500 T T T 0.5000 0.2500 0.7500 T T T 0.7500 0.5000 0.2500 T T T 0.7500 0.2500 0.7500 T T T FIN echo "a=$a ncpu=$ncpu" cmd="mpiexec -np $ncpu vasp" $cmd >& vasp.log E=`tail -1 OSZICAR` echo $a $E | sed -s 's/F=//; s/E0=//; s/d E =//;' >> Elatt.C44.dat done
THIS SCRIPT IS NOT A COMPLETED ONE
This script calculates C44 employing the shear method
Script for C44 shear method
#!/bin/bash rm WAVECAR for a in -0.003 -0.002 -0.001 0.00 0.001 0.002 0.003 do cat > POSCAR << FIN POSCAR for FCC Si 5.402856 universal scaling factor 1.0 0.0 0.0 first Bravais lattice vector $a 1.0 0.0 second Bravais lattice vector 0.0 0.0 1.0 third Bravais lattice vector 8 number of atoms per species selective dynamics direct direct or cart (only first letter is significant) 0.00000 0.00000 0.00000 T T T 0.50000 0.50000 0.00000 T T T 0.00000 0.50000 0.50000 T T T 0.50000 0.00000 0.50000 T T T 0.25000 0.25000 0.25000 T T T 0.75000 0.75000 0.25000 T T T 0.25000 0.75000 0.75000 T T T 0.75000 0.25000 0.75000 T T T FIN
Analyze data
After running the two scripts above for C11 and C44, the data files: Elatt.C11.dat and Elatt.C44.dat will be created.
If you choose your POTCAR for US-LDA PP, the date should be similar to
Elatt.C11.dat
0.997 1 -.47813207E+02 -.47813293E+02 -.478132E+02 0.998 1 -.47813604E+02 -.47813690E+02 -.478136E+02 0.999 1 -.47813844E+02 -.47813929E+02 -.478138E+02 1.000 1 -.47813926E+02 -.47814012E+02 -.478139E+02 1.001 1 -.47813851E+02 -.47813936E+02 -.478139E+02 1.002 1 -.47813619E+02 -.47813705E+02 -.478136E+02 1.003 1 -.47813231E+02 -.47813316E+02 -.478132E+02
Elatt.C44.dat
0.997 3 -.95626284E+02 -.95626386E+02 -.264643E-03 0.998 3 -.95627212E+02 -.95627315E+02 -.116818E-03 0.999 3 -.95627770E+02 -.95627874E+02 -.289956E-04 1.000 1 -.95627959E+02 -.95628064E+02 -.956280E+02 1.001 3 -.95627782E+02 -.95627888E+02 -.285597E-04 1.002 3 -.95627238E+02 -.95627345E+02 -.113344E-03 1.003 3 -.95626329E+02 -.95626438E+02 -.252950E-03
Download the following Matlab functions. fit_a0EB.m, fit_C11EB.m, fit_C44EB.m,
Place the data files of Elatt.C11.dat, Elatt.C44.dat, Elatt.B.dat from Bulk modulus calculation of Si
Launch Matlab, and type the following commands.
[a0, Ecoh, B] = fit_a0EB('Elatt.B.dat',2,8);
[C11,C12] = fit_C11EB('Elatt.C11.dat',a0,B);
C44 = fit_C44EB('Elatt.C44.dat',a0,C11,C12);
The first line fits the energy data to a quadratic curve and computes the equilibrium lattice constant, cohesive energy and bulk modulus([a0, Ecoh, B]) respectively. The parameter 2 specifies the number of atoms in the computation cell and 8 is the number of atoms in a conventional unit cell. The second lines computes the C11 and C12 elastic constants and the third line calculates C44. Suggested results are C11 = 161 GPa, C12 = 63 GPa and C44 = 76 GPa for calculations with US-LDA PP.
C44 data we got here is much smaller than the paper's result of 105 GPa but is the same as the experimental value of 76 GPa. We suspect maybe 105 GPa is the unrelaxed result since only specifying ISIF=2 doesn't enable relaxation.