calculate the dipole moment for rippled GeSe2 structure.
├── 0.02500_w_H
│ └── scf
│ ├── cfit.py
│ ├── OUTCAR
│ └── POSCAR
├── 0.02500_wo_H
│ └── scf
│ ├── OUTCAR
│ └── POSCAR
├── 0.05000_w_H
│ └── scf
│ ├── cfit.py
│ ├── OUTCAR
│ └── POSCAR
├── 0.05000_wo_H
│ └── scf
│ ├── OUTCAR
│ └── POSCAR
├── 0.07500_w_H
│ └── scf
│ ├── cfit.py
│ ├── OUTCAR
│ └── POSCAR
├── 0.07500_wo_H
│ └── scf
│ ├── OUTCAR
│ └── POSCAR
├── 0.10000_w_H
│ └── scf
│ ├── cfit.py
│ ├── OUTCAR
│ └── POSCAR
├── 0.10000_wo_H
│ └── scf
│ ├── OUTCAR
│ └── POSCAR
├── cfit.py
├── p.dat
├── process.ipynb
├── proc.py
├─md.py
├─feature.json
└── proc.sh
There are ten folders that used for VASP dipole moment calculation. The name of every folder consist three parts: strain_flat/rippled_Hydrogen, for example:
the 0.10000_wo_H
folder stands for flat GeSe2 structure with 0.1 strain and saturated by Hydrogen.
the 0.10000_w_H
folder stands for rippled GeSe2 structure with 0.1 strain and saturated by Hydrogen.
The circle_fit module supply two methods to calcute the curvature radius. In our work, we use the average value.
from pymatgen import Structure
#from circle_fit import hyper_fit, least_squares_circle
import circle_fit as cf
import copy
st=Structure.from_file('POSCAR')
sites=[site for site in st if site.species_string=='Ge']
sites=sites[1:-1]
X=[site.coords[0] for site in sites]
Y=[site.coords[2] for site in sites]
x=[]
y=[]
for i in range(len(X)):
x.append([X[i]])
y.append([Y[i]])
coords = [[x[i][0], y[i][0]] for i in range(len(x))]
xc0,yc1,r0,s0 = cf.hyper_fit(coords)
xc1,yc1,r1,s1 = cf.least_squares_circle(coords)
print((r0+r1)/2)
path=`pwd`
list='0.02500 0.05000 0.07500 0.10000'
rm -f p.dat
for i in $list
do
cd $path/${i}_wo_H/scf
ion1=`grep 'p\[ion\]=' OUTCAR | cut -d '=' -f 2 | awk '{print $2,$3,$4}'`
elc1=`grep 'p\[elc\]=' OUTCAR | cut -d '=' -f 2 | awk '{print $2,$3,$4}'`
latc1=`sed -n 5p POSCAR | awk '{print $3}'`
e0=`grep TOTEN OUTCAR | tail -1 | awk '{ print $5}'`
cd $path/${i}_w_H/scf
ion2=`grep 'p\[ion\]=' OUTCAR | cut -d '=' -f 2 | awk '{print $2,$3,$4}'`
elc2=`grep 'p\[elc\]=' OUTCAR | cut -d '=' -f 2 | awk '{print $2,$3,$4}'`
latc2=`sed -n 5p POSCAR | awk '{print $3}'`
cp $path/cfit.py .
rc=`python cfit.py`
e1=`grep TOTEN OUTCAR | tail -1 | awk '{print $5}'`
echo $i $ion1 $elc1 $ion2 $elc2 $latc1 $latc2 $rc $e0 $e1>>$path/p.dat
cd $path
done
cat >proc.py <<EOF
import numpy as np
ret=np.loadtxt('p.dat')
for i in range(ret.shape[0]):
ptwo=ret[i,3]+ret[i,6]
latc1=ret[i,-5]
ptw=ret[i,9]+ret[i,12]
latc2=ret[i,-4]
r=ret[i,0]
de=ret[i,-1]-ret[i,-2]
print("strain %.3f dipole_flat %.3f dipole_rippled %.3e r %.3f e %.3f"%(r,ptwo/latc1,ptw%latc2*1.602e-19/(ret[i,-3]*1e-10)/1e-12,ret[i,-3],de))
EOF
python proc.py
After running sh proc.sh
, we can directly obtain the results:
strain 0.025 dipole_flat -0.000 dipole_rippled 3.311e-01 r 165.767 e 0.045
strain 0.050 dipole_flat -0.000 dipole_rippled 9.065e-01 r 94.109 e 0.026
strain 0.075 dipole_flat -0.000 dipole_rippled 1.575e+00 r 71.903 e -0.030
strain 0.100 dipole_flat -0.000 dipole_rippled 1.918e+00 r 63.701 e -0.111
where, r
stands for curvature radius and e
stands for energy difference between flat structure and rippled one.
install the pymatgen and deepchem package
pip install pymatgen
pip install deepchem
run the script
python md.py
此处可能存在不合适展示的内容,页面不予展示。您可通过相关编辑功能自查并修改。
如您确认内容无涉及 不当用语 / 纯广告导流 / 暴力 / 低俗色情 / 侵权 / 盗版 / 虚假 / 无价值内容或违法国家有关法律法规的内容,可点击提交进行申诉,我们将尽快为您处理。