VASP Computing Elastic Constants of Si: Difference between revisions

From Micro and Nano Mechanics Group
Jump to navigation Jump to search
 
(27 intermediate revisions by 2 users not shown)
Line 3: Line 3:
VASP: Computing Elastic Constants of Si</STRONG></font></P>
VASP: Computing Elastic Constants of Si</STRONG></font></P>
<DIV>
<DIV>
<P ALIGN="CENTER"><STRONG> </STRONG></P>
<P ALIGN="CENTER"><STRONG>Adriano Sanchez，Yanming Wang and Wei Cai </STRONG></P>
</DIV>
</DIV>


Line 9: Line 9:
===Input files===
===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 <tt>MC2</tt> in PARALLEL mode in <tt>~/Codes/VASP/runs/Si_bulk_C11</tt> and <tt>~/Codes/VASP/runs/Si_bulk_C44</tt> two directories. Each directory contains the following files.
'''(Under construction!!!)'''

Here we give an example of how to use VASP to compute the elastic constants C11 and C44 of LDA-Si. We performed this calculation on <tt>MC-CC</tt> in PARALLEL model in the <tt>~/Codes/VASP/runs/Au/LDA/perfect.21x21x21</tt> directory. This directory contains the following files.


'''INCAR'''
'''INCAR'''


<pre>
<pre>
PREC = High
ENCUT = 400 or 500??

ISIF = 2 It relaxes atoms but keep cell volume fixed!
ISMEAR = 1
ISTART = 0
ICHARG = 2
ISMEAR = 1
SIGMA = 0.1
SIGMA = 0.1

EDIFF = 1E-09
NELM = 40
ENMAX = 500
ENCUT = 500

ISIF = 2
NSW = 100
IBRION = 2
</pre>
</pre>


Line 32: Line 41:
</pre>
</pre>


'''POSCAR example for C11'''
'''POSCAR '''


Here's an example for Si C11 and C12 calculation, which will be changed later in PBS script every iteration.
<pre>
<pre>
POSCAR for DC Si
POSCAR for DC Si
Line 53: Line 63:
</pre>
</pre>


Here's an example for Si C44 calculation, which will be changed later in PBS script every iteration.
To do this calculation, you also need to put the LDA pseudopotential file as <tt>POTCAR</tt> in this directory.

'''POSCAR example for C44'''


<pre>
<pre>
Line 83: Line 91:
1.0000 0.2500 0.7500 T T T
1.0000 0.2500 0.7500 T T T
</pre>
</pre>

'''POTCAR'''

In this example, we choose the ultrasoft-LDA PP of Si provided by VASP.


===Run VASP===
===Run VASP===


To compute the elastic constants we use the following script <tt>auto.B.serial</tt> to run <tt>vasp</tt> repeated with different values for the first C11 component of the first Bravais lattice vector, '''$a'''. Make sure you have previously calculated the equilibrium lattice constant for this supercell with 4 atoms, otherwise errors in the calculations could arise. The script creates automatically the POSCAR file as in the example above for every '''a''' especified in the script in the range from 0.997 to 1.003
To compute the elastic constants, we prepared the following PBS scripts to submit the jobs to MC2 cluster. <tt>vasp.pbs.C11</tt> is used to to run <tt>vasp</tt> 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 <tt>vasp.pbs.C44</tt> changes the value of '''$a''' for the first and second component of the first Bravais lattice vector.


'''Script for C11'''
'''PBS script for C11'''


<pre>
<pre>
#!/bin/bash
#!/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
rm WAVECAR
for a in 0.997 0.998 0.999 1.000 1.001 1.002 1.003
for a in 0.997 0.998 0.999 1.000 1.001 1.002 1.003
Line 113: Line 152:
0.7500 0.2500 0.7500 T T T
0.7500 0.2500 0.7500 T T T
FIN
FIN

echo "a=$a"
echo "a=$a ncpu=$ncpu"
./vasp

cmd="mpiexec -np $ncpu vasp"
$cmd >& vasp.log


E=`tail -1 OSZICAR`
E=`tail -1 OSZICAR`
echo $a $E | sed -s 's/F=//; s/E0=//; s/d E =//;' >> Elatt.C11.dat
echo $a $E | sed -s 's/F=//; s/E0=//; s/d E =//;' >> Elatt.C11.dat

p=`grep pressure OUTCAR | cut -b 25-34`
echo $a $p >> platt.C11.dat


done
done
</pre>
</pre>


'''Script for C44'''
'''PBS script for C44'''


<pre>
<pre>
#!/bin/bash
#!/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
rm WAVECAR
for a in 0.997 0.998 0.999 1.000 1.001 1.002 1.003
for a in 0.997 0.998 0.999 1.000 1.001 1.002 1.003
Line 150: Line 216:
0.7500 0.2500 0.5000 T T T
0.7500 0.2500 0.5000 T T T
0.2500 0.5000 0.2500 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.5000 0.7500 0.7500 T T T
0.2500 0.0000 0.2500 T T T
0.2500 0.0000 0.2500 T T T
Line 155: Line 223:
0.5000 0.2500 0.7500 T T T
0.5000 0.2500 0.7500 T T T
0.7500 0.5000 0.2500 T T T
0.7500 0.5000 0.2500 T T T
1.0000 0.7500 0.7500 T T T
0.7500 0.2500 0.7500 T T T
1.0000 0.2500 0.7500 T T T
FIN
FIN

echo "a=$a"
echo "a=$a ncpu=$ncpu"
./vasp

cmd="mpiexec -np $ncpu vasp"
$cmd >& vasp.log


E=`tail -1 OSZICAR`
E=`tail -1 OSZICAR`
echo $a $E | sed -s 's/F=//; s/E0=//; s/d E =//;' >> Elatt.C44.dat
echo $a $E | sed -s 's/F=//; s/E0=//; s/d E =//;' >> Elatt.C44.dat

p=`grep pressure OUTCAR | cut -b 25-34`
echo $a $p >> platt.C44.dat


done
done
</pre>

'''THIS SCRIPT IS NOT A COMPLETED ONE'''

This script calculates C44 employing the shear method

'''Script for C44 shear method'''

<pre>
#!/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
</pre>
</pre>


===Analyze data===
===Analyze data===


After running the two scripts above for C11 and C44, <tt>./auto.B.serial</tt> it will create the data files: <tt>Elatt.C11.dat</tt>,<tt>platt.C11.dat</tt>, <tt>Elatt.C44.dat</tt>, and <tt>platt.C44.dat</tt>.
After running the two scripts above for C11 and C44, the data files: <tt>Elatt.C11.dat</tt> and <tt>Elatt.C44.dat</tt> will be created.


If you choose your POTCAR for US-LDA PP, the date should be similar to


<tt>Elatt.C11.dat</tt>
Launch <tt>octave</tt> and run the following functions [[media:fit_a0EB.m.txt | fit_a0EB.m]] and [[media:fit_a0B.m.txt | fit_a0B.m]], '''cal_VASP.m'''
<pre>
fit_a0EB('Elatt.B.dat',4,4);
0.997 1 -.47813207E+02 -.47813293E+02 -.478132E+02
fit_C11EB('Elatt.C11.dat',a0,B);
0.998 1 -.47813604E+02 -.47813690E+02 -.478136E+02
fit_C44EB('Elatt.C44.dat',a0,C11,C12);
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
</pre>


<tt>Elatt.C44.dat</tt>
In MATLAB click "change folder" when you open the file and remove the semicolon to run it from the command line.
<pre>
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.
0.997 3 -.95626284E+02 -.95626386E+02 -.264643E-03
Refer to previous calculation?!)
0.998 3 -.95627212E+02 -.95627315E+02 -.116818E-03
The second lines computes the C11 and C12 elastic constants and the third line calculates C44. You should obtain C11 = 161 GPa, C12 = 179.2 GPa and C44 = 42 GPa.
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
</pre>


Download the following Matlab functions.
===Parallel computation===
[[media:fit_a0EB.m.txt | fit_a0EB.m]],
[[media:fit_C11EB.m.txt | fit_C11EB.m]],
[[media:fit_C44EB.m.txt | fit_C44EB.m]],


Place the data files of <tt>Elatt.C11.dat</tt>, <tt>Elatt.C44.dat</tt>, <tt>Elatt.B.dat</tt> from [[VASP_Computing_Bulk_Modulus_of_Si | Bulk modulus calculation of Si]]
To run <tt>vasp</tt> in parallel, you need to submit [[media:vasp.pbs.txt | vasp.pbs]] as


Launch <tt>Matlab</tt>, and type the following commands.
qsub vasp.pbs
[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.


<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>
You will need the following two files to do this calculation in parallel on SU-AHPCRC:
[[media:auto.B.par.Au.txt | auto.B.par]] and
[[media:B.pbs.Au.txt | B.pbs]].

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.

Retrieved from "http://micro.stanford.edu/mediawiki/index.php?title=VASP_Computing_Elastic_Constants_of_Si&oldid=6309"