2D材料方向依赖的泊松比和杨氏模量图

materials

发布日期: 2021-07-01

文章字数: 1.1k

阅读时长: 6 分

阅读次数:

问题描述

2D材料的杨氏模量和泊松比是材料非常重用的力学属性，我们如何通过vasp的计算来画出其方向依赖的泊松比和杨氏模量图：

数据提取

关于弹性常数的计算可以采用如下的INCAR参数，通过vasp直接计算（也可以通过手动应力-应变的方式来计算）

# INCAR
SYSTEM = system
PREC   = A
ISPIN  =  1
LREAL  = Auto
ALGO  =  Normal

#Electronic Relaxation
ENCUT  =  500
NELM   =  200
NELMIN=  6
NELMDL= -5
EDIFF  = 1E-6
#LVDW = TRUE

#Ionic relaxation
EDIFFG = -0.005
NSW    =  1
IBRION =  6
ISIF   =  3
ISYM   =  2
NFREE  =  4
POTIM  =  0.01

#DOS related values
ISMEAR = 0
SIGMA  = 0.05

#Write flags
LWAVE  =  F
LCHARG =  F

计算完成之后，在当前文件夹下执行如下的bash脚本(sh mesh.sh 3.0)，完成数据的提取和处理

注：脚本默认c方向为真空层方向，建模时需要留意；此外，3.0的数值表示材料的实际有效厚度，需要根据自己的材料测量设置。

c11=`grep -A 9 "TOTAL ELASTIC" OUTCAR | tail -n 8 | grep XX | awk '{print $2}'`
c22=`grep -A 9 "TOTAL ELASTIC" OUTCAR | tail -n 8 | grep YY | awk '{print $3}'`
c12=`grep -A 9 "TOTAL ELASTIC" OUTCAR | tail -n 8 | grep XX | awk '{print $3}'`
c66=`grep -A 9 "TOTAL ELASTIC" OUTCAR | tail -n 8 | grep XY | awk '{print $5}'`

echo "-----------Elastic Constants are-------------"
C11=`echo "scale=2;$c11/10.0" | bc`
C22=`echo "scale=2;$c22/10.0" | bc`
C12=`echo "scale=2;$c12/10.0" | bc`
C66=`echo "scale=2;$c66/10.0" | bc`
echo "C11 :  $C11 "
echo "C22 :  $C22 "
echo "C12 :  $C12 "
echo "C66 :  $C66 "
echo "---------------------------------------------"

res=`echo "$C11 >0" | bc`
if [ $res == 1 ];then
   echo "C11 > 0"
else
   echo "C11 < 0"
fi

res=`echo "$C22 >0" | bc`
if [ $res == 1 ];then
   echo "C22 > 0"
else
   echo "C22 < 0"
fi

res=`echo "$C66 >0" | bc`
if [ $res == 1 ];then
   echo "C66 > 0"
else
   echo "C66 < 0"
fi

res=`echo "$C11*$C22-$C12*$C12 >0" | bc`
if [ $res == 1 ];then
   echo "C11xC22-C12**2 > 0"
else
   echo "C11xC22-C12**2 < 0"
fi

c=`sed -n 5p POSCAR | awk '{print $3}'`
cat>proc.py<<!
import numpy as np
from numpy import sin,cos,pi
L0=$1;
c=$c;
c_L0=c/L0;
C11=$C11*c_L0;
C22=$C22*c_L0;
C12=$C12*c_L0;
C66=$C66*c_L0;
theta=np.linspace(0,2*pi,360);
v_zz=C12/C22;
d1=C11/C22+1-(C11*C22-C12**2)/C22/C66;
d2=-(2*C12/C22-(C11*C22-C12**2)/C22/C66);
d3=C11/C22;
Y_zz=(C11*C22-C12**2)/C22;
E=Y_zz/((cos(theta))**4+d2*(cos(theta))**2.*(sin(theta))**2+d3*(sin(theta))**4);
V=(v_zz*(cos(theta))**4-d1*(cos(theta))**2.*(sin(theta))**2+v_zz*(sin(theta))**4)/((cos(theta))**4+d2*(cos(theta))**2.*(sin(theta))**2+d3*(sin(theta))**4);
s=sin(theta);
c=cos(theta);
E1=1./((C22*c**4-2*C12*c**2.*s**2+C11*s**4)/(C11*C22-C12**2)+c**2.*s**2/C66);
V1=(C12*c**4-(C11+C22-(C11*C22-C12**2)/C66)*c**2.*s**2+C12*s**4)/(C22*c**4-(2*C12-(C11*C22-C12**2)/C66)*c**2.*s**2+C11*s**4);
res=np.vstack((theta,E,V)).T
np.savetxt('EV.dat',res,fmt='%.3f %.3f %.3f')
print("%.3f %.3f %.3f %.3f"%(C11,C12,C22,C66))
!
# call the python to calculate the Young's modulus and Poisson's ratio.
python < proc.py > log.dat

可以通过log.dat文件来查看是否满足力学稳定性条件以及C11,C12,C22,C66的值。

画图

在当前文件夹下保存mech_polar.py

import matplotlib.pyplot as plt
import numpy as np
fn="EV.dat"
f=np.loadtxt(fn)

fig = plt.figure()

ax=fig.add_subplot(121, projection='polar')# ,facecolor="lightgoldenrodyellow")
skip=2
theta= f[0::skip,0]
Y=f[0::skip,1]
ax.plot(theta, Y, color="tab:orange", lw=1, ls="--", marker='h',alpha=0.6, label="$E$")
ax.set_rlabel_position(10)  # get radial labels away from plotted line
ax.tick_params(grid_color="palegoldenrod")
angle = np.deg2rad(67.5)
ax.legend(loc="lower left", bbox_to_anchor=(.5 + np.cos(angle)/2, .5 + np.sin(angle)/2))

ax1=fig.add_subplot(122, projection='polar') #,facecolor="lightgoldenrodyellow")
ax1.tick_params(grid_color="palegoldenrod")
V=f[:,2]
T=f[:,0]
Vp_idx=np.where(V>0)
Vn_idx=np.where(V<0)

skip1=3
Vp=V[Vp_idx][0::skip1]
theta_p=T[Vp_idx][0::skip1]

skip2=1
Vn=V[Vn_idx][0::skip2]
theta_n=T[Vn_idx][0::skip2]

ax1.plot(theta_p, Vp, color="tab:green", lw=1, ls="--", marker='h',alpha=0.6, label=r"+ $\nu$")
ax1.plot(theta_n, Vn, color="tab:red", lw=0.1, marker='o',alpha=.5, label=r"- $\nu$")
angle = np.deg2rad(67.5)
ax1.legend(loc="lower left", bbox_to_anchor=(.5 + np.cos(angle)/2, .5 + np.sin(angle)/2))
#ax1.set_rmax(2)
ax1.set_rlabel_position(45)  # get radial labels away from plotted line
ax1.set_rticks([0.3, 0.60, 0.9])  # less radial ticks

plt.tight_layout()
#plt.savefig('EV.jpg',format='jpg',dpi=300)
plt.show()