VASP Computing Ideal Shear Strength of Au: Difference between revisions
| Line 193: | Line 193: | ||
proc run_vasp { ncpu } { |
proc run_vasp { ncpu } { |
||
set fileID [open "run_vasp" w ] |
|||
#if { $ncpu == 1 } { |
|||
| ⚫ | |||
| ⚫ | |||
| ⚫ | |||
#} else { |
|||
file attributes run_vasp -permissions u+x |
|||
| ⚫ | |||
| ⚫ | |||
| ⚫ | |||
file attributes run_vasp -permissions u+x |
|||
exec ./run_vasp |
|||
#} |
|||
} |
} |
||
Revision as of 06:55, 15 February 2009
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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \gamma_{xy}} in our coordinate system. Due to elastic anisotropy, we need to adjust other strain components, so that the only non-zero stress component is Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sigma_{xy}} . Due to the symmetry of this problem, it turns out that we only need to adjust normal strain components.
For each applied strain Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \gamma_{xy}} , we need to run a series of VASP calculations to find the strain components that relax all other stress components except Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sigma_{xy}} . 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 } {
set fileID [open "run_vasp" w ]
puts $fileID "/opt/mpiexec/bin/mpiexec --comm=pmi -np $ncpu ../../../../bin/vasp.caiwei.mva2"
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