VASP Computing Ideal Shear Strength of Au
Jump to navigation
Jump to search
VASP: Ideal Shear Strength of Au
Input files
Here are the basic input files required for VASP calculation. Some of the files need to be changed since we need to perform a large number of calculations.
INCAR
System = fcc Au LWAVE = .FALSE. ENCUT = 400 ISMEAR = 1 SIGMA = 0.2 ISIF = 2
To increase plane wave cutoff, we manually put PREC = HIGH or PREC = Accurate. Default value is PREC = Medium
KPOINTS
43x43x43 0 0 = automatic generation of k-points Monkhorst 43 43 43 0 0 0
POSCAR
POSCAR for FCC Au (created by tcl) 4.0615 0.612372435696 0.000000000000 -0.353553390593 0.408248290464 0.577350269190 0.000000000000 0.612372435696 0.000000000000 0.353553390593 1 Cartesian (real coordinates r) 0 0 0
You also need to put the LDA pseudopotential file as POTCAR in this directory.
Results
Perfect crystal
First, we perform a test simulation in which we relax the electrons in a gold perfect crystal.
qsub vasp.pbs
The vasp.pbs script is listed below.
#!/bin/bash #PBS -N Au_perf.43.E4.8 #PBS -j oe #PBS -l nodes=1:ppn=8,walltime=24: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 /opt/mpiexec/bin/mpiexec --comm=pmi -np $ncpu ../../../../bin/vasp.caiwei.mva2 >& $PBS_O_WORKDIR/vasp.log
We run this test case on SU-AHPCRC. The calculation takes 160 seconds (on 8 CPUs) to finish. The ground state energy is E = -4.3951 eV.
Shear deformation
Next, we apply shear strain on (111) plane along the [1 -1 0] direction. This coresponds to in our coordinate system. Due to elastic anisotropy, we need to adjust other strain components, so that the only non-zero stress component is . Due to the symmetry of this problem, it turns out that we only need to adjust normal strain components.
For each applied strain , we need to run a series of VASP calculations to find the strain components that relax all other stress components except . After that, is increased a little further and the relaxation is repeated. For each , the shear stress at the end of the relaxation is recorded. The ideal shear stress is the maximum recorded value for .
The entire calculation is submitted by the following command.
qsub tcl.pbs
The tcl.pbs script is listed below.
#!/bin/bash #PBS -N Au_shear.43.E4.rlx #PBS -j oe #PBS -l nodes=1:ppn=8,walltime=96: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 /usr/bin/tclsh shear_rlx.tcl $ncpu
The shear_rlx.tcl script is listed below.
#####################################
# definition of functions
#####################################
proc write_POSCAR { posfile a0 ax ay az bx by bz cx cy cz } {
#set posfile "POSCAR"
set fileID [open $posfile w ]
puts $fileID "POSCAR for FCC Au (created by tcl)"
puts $fileID "$a0"
puts $fileID [format "%15.12f %15.12f %15.12f" $ax $ay $az]
puts $fileID [format "%15.12f %15.12f %15.12f" $bx $by $bz]
puts $fileID [format "%15.12f %15.12f %15.12f" $cx $cy $cz]
puts $fileID "1"
puts $fileID "Cartesian (real coordinates r)"
puts $fileID "0 0 0"
close $fileID
}
proc read_POSCAR { posfile } {
#set posfile "POSCAR"
set fileID [open $posfile r ]
# first line is comment
gets $fileID line
# second line is a0
gets $fileID line
scan $line "%f" a0
# now comes three lines specifying the cell
gets $fileID line
scan $line "%f %f %f" ax ay az
gets $fileID line
scan $line "%f %f %f" bx by bz
gets $fileID line
scan $line "%f %f %f" cx cy cz
# skip the rest of the file
close $fileID
# return cell geometry as a list
return [list $a0 $ax $ay $az $bx $by $bz $cx $cy $cz]
}
proc get_stress { } {
set fileID [open "get_stress" w ]
puts $fileID "grep \"in kB\" OUTCAR | sed \'s/in kB //\' "
close $fileID
file attributes get_stress -permissions u+x
exec ./get_stress
}
proc run_vasp { ncpu } {
#if { $ncpu == 1 } {
# exec ./vasp
#} else {
set fileID [open "run_vasp" w ]
puts $fileID "mpiexec --comm=ib -np $ncpu ../../../../bin/vasp.2008-11-15"
close $fileID
file attributes run_vasp -permissions u+x
exec ./run_vasp
#}
}
#####################################
# beginning of main program
#####################################
if { $argc == 0 } {
set ncpu 1
} elseif { $argc > 0 } {
set ncpu [lindex $argv 0]
}
puts "ncpu = $ncpu"
# modulus in kilobar (100 MPa)
set C11 2020
set C12 1690
set C44 450
if { 0 } {
# compute nu, E using isotropic approximation
set lambda $C12
set mu [expr ($C11-$C12)/2.0 ]
set nu [expr $lambda / (2.0*($lambda + $mu)) ]
# lambda = 2*mu*nu / (1-2*nu)
set YoungsModulus [expr 2*$mu*(1+$nu) ]
puts "mu = $nu nu = $nu E = $YoungsModulus"
} else {
# directly set nu, E
set YoungsModulus 780
set nu 0.44
puts "nu = $nu E = $YoungsModulus"
}
set S11 [expr 1/$YoungsModulus ]
set S12 [expr -$nu/$YoungsModulus ]
puts "S11 = $S11 S12 = $S12"
set factor 0.4
set sig_tol 0.1
set a0 4.0615
#set a0 4.062
set ax0 [expr sqrt(6)/4]
set ay0 0
set az0 [expr -sqrt(2)/4]
set bx0 [expr sqrt(6)/6]
set by0 [expr sqrt(3)/3]
set bz0 0
set cx0 [expr sqrt(6)/4]
set cy0 0
set cz0 [expr sqrt(2)/4]
write_POSCAR "POSCAR.0" $a0 $ax0 $ay0 $az0 $bx0 $by0 $bz0 $cx0 $cy0 $cz0
#run_vasp $ncpu
puts "bx0 = $bx0"
set datafile "stress_rlx.dat"
set rlxfileID [open $datafile w ]
set max_dx_iter 50
set stop_dx_iter 50
set start_dx_iter 1
set reset_poscar 1
set max_rlx_iter 20
for { set dx_iter $start_dx_iter } { $dx_iter <= $stop_dx_iter } { incr dx_iter 1 } {
# set cell dimension: bx
set dx [expr sqrt(6)/6 * $dx_iter/$max_dx_iter]
set bx [expr $bx0 + $dx]
puts "dx = $dx bx = $bx"
# set cell dimensions other than bx
if { $dx_iter == $start_dx_iter } {
if { $reset_poscar } {
# initialize cell value
set ax $ax0
set cx $cx0
set by $by0
set az $az0
set cz $cz0
write_POSCAR "POSCAR" $a0 $ax $ay0 $az $bx $by $bz0 $cx $cy0 $cz
write_POSCAR "POSCAR.$dx_iter" $a0 $ax $ay0 $az $bx $by $bz0 $cx $cy0 $cz
} else {
# read cell value from POSCAR
set cell [ read_POSCAR "POSCAR" ]
set a0 [lindex $cell 0 ]
set ax [lindex $cell 1 ]
set ay [lindex $cell 2 ]
set az [lindex $cell 3 ]
# set bx [lindex $cell 4 ]
set by [lindex $cell 5 ]
set bz [lindex $cell 6 ]
set cx [lindex $cell 7 ]
set cy [lindex $cell 8 ]
set cz [lindex $cell 9 ]
puts "read cell geometry:"
puts " $a0"
puts " $ax $ay $az"
puts " $bx $by $bz"
puts " $cx $cy $cz"
}
}
# adjusting box to relax normal stress components to zero
for { set rlx_iter 0 } { $rlx_iter <= $max_rlx_iter } { incr rlx_iter 1 } {
run_vasp $ncpu
set stress [ get_stress ]
scan $stress "%f %f %f %f %f %f" sig_xx sig_yy sig_zz sig_xy sig_yz sig_zx
puts $rlxfileID "$dx $rlx_iter $stress"
flush $rlxfileID
if { $rlx_iter < $max_rlx_iter } {
if { 0 } {
# account for Poisson's effect
set e_xx [expr $sig_xx * $S11 + $sig_yy * $S12 + $sig_zz * $S12]
set e_yy [expr $sig_xx * $S12 + $sig_yy * $S11 + $sig_zz * $S12]
set e_zz [expr $sig_xx * $S12 + $sig_yy * $S12 + $sig_zz * $S11]
set bx [expr $e_yy * ($bx0+$dx) * $factor + $bx ]
} else {
# do not account for Poisson's effect
set e_xx [expr $sig_xx / $C11 ]
set e_yy [expr $sig_yy / $C11 ]
set e_zz [expr $sig_zz / $C11 ]
}
puts "sig_xx = $sig_xx sig_yy = $sig_yy sig_zz = $sig_zz"
puts "e_xx = $e_xx e_yy = $e_yy e_zz = $e_zz"
set ax [expr $e_xx * $ax0 * $factor + $ax ]
set cx [expr $e_xx * $cx0 * $factor + $cx ]
set by [expr $e_yy * $by0 * $factor + $by ]
set az [expr $e_zz * $az0 * $factor + $az ]
set cz [expr $e_zz * $cz0 * $factor + $cz ]
write_POSCAR "POSCAR" $a0 $ax $ay0 $az $bx $by $bz0 $cx $cy0 $cz
write_POSCAR "POSCAR.$dx_iter.[expr $rlx_iter+1]" $a0 $ax $ay0 $az $bx $by $bz0 $cx $cy0 $cz
}
}
}
close $rlxfileID
# End of program
The TCL script uses the run_vasp script, which is listed below.
/opt/mpiexec/bin/mpiexec --comm=pmi -np 8 ../../../../bin/vasp.caiwei.mva2