Skip to content

Commit 7c4ba3b

Browse files
MabinogiXMabinogiX
committed
代码重构
1 parent fd45ceb commit 7c4ba3b

23 files changed

+576
-418
lines changed

‎IR-intensity.sh

+105
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
#!/bin/bash
2+
# A utility for calculating the vibrational intensities from VASP output (OUTCAR)
3+
# (C) David Karhanek, 2011-03-25, ICIQ Tarragona, Spain (www.iciq.es)
4+
5+
# extract Born effective charges tensors
6+
printf "..reading OUTCAR"
7+
BORN_NROWS=`grep NIONS OUTCAR | awk '{print $12*4+1}'`
8+
if [ `grep 'BORN' OUTCAR | wc -l` = 0 ] ; then \
9+
printf " .. FAILED! Born effective charges missing! Bye! \n\n" ; exit 1 ; fi
10+
grep "in e, cummulative" -A $BORN_NROWS OUTCAR > born.txt
11+
12+
# extract Eigenvectors and eigenvalues
13+
if [ `grep 'SQRT(mass)' OUTCAR | wc -l` != 1 ] ; then \
14+
printf " .. FAILED! Restart VASP with NWRITE=3! Bye! \n\n" ; exit 1 ; fi
15+
EIG_NVIBS=`grep -A 2000 'SQRT(mass)' OUTCAR | grep 'cm-1' | wc -l`
16+
EIG_NIONS=`grep NIONS OUTCAR | awk '{print $12}'`
17+
EIG_NROWS=`echo "($EIG_NIONS+3)*$EIG_NVIBS+3" | bc`
18+
grep -A $(($EIG_NROWS+2)) 'SQRT(mass)' OUTCAR | tail -n $(($EIG_NROWS+1)) | sed 's/f\/i/fi /g' > eigenvectors.txt
19+
printf " ..done\n"
20+
21+
# set up a new directory, split files - prepare for parsing
22+
printf "..splitting files"
23+
mkdir intensities ; mv born.txt eigenvectors.txt intensities/
24+
cd intensities/
25+
let NBORN_NROWS=BORN_NROWS-1
26+
let NEIG_NROWS=EIG_NROWS-3
27+
let NBORN_STEP=4
28+
let NEIG_STEP=EIG_NIONS+3
29+
tail -n $NBORN_NROWS born.txt > temp.born.txt
30+
tail -n $NEIG_NROWS eigenvectors.txt > temp.eige.txt
31+
mkdir inputs ; mv born.txt eigenvectors.txt inputs/
32+
split -a 3 -d -l $NEIG_STEP temp.eige.txt temp.ei.
33+
split -a 3 -d -l $NBORN_STEP temp.born.txt temp.bo.
34+
mkdir temps01 ; mv temp.born.txt temp.eige.txt temps01/
35+
for nu in `seq 1 $EIG_NVIBS` ; do
36+
let nud=nu-1 ; ei=`printf "%03u" $nu` ; eid=`printf "%03u" $nud` ; mv temp.ei.$eid eigens.vib.$ei
37+
done
38+
for s in `seq 1 $EIG_NIONS` ; do
39+
let sd=s-1 ; bo=`printf "%03u" $s` ; bod=`printf "%03u" $sd` ; mv temp.bo.$bod borncs.$bo
40+
done
41+
printf " ..done\n"
42+
43+
# parse deviation vectors (eig)
44+
printf "..parsing eigenvectors"
45+
let sad=$EIG_NIONS+1
46+
for nu in `seq 1 $EIG_NVIBS` ; do
47+
nuu=`printf "%03u" $nu`
48+
tail -n $sad eigens.vib.$nuu | head -n $EIG_NIONS | awk '{print $4,$5,$6}' > e.vib.$nuu.allions
49+
split -a 3 -d -l 1 e.vib.$nuu.allions temp.e.vib.$nuu.ion.
50+
for s in `seq 1 $EIG_NIONS` ; do
51+
let sd=s-1; bo=`printf "%03u" $s`; bod=`printf "%03u" $sd`; mv temp.e.vib.$nuu.ion.$bod e.vib.$nuu.ion.$bo
52+
done
53+
done
54+
printf " ..done\n"
55+
56+
# parse born effective charge matrices (born)
57+
printf "..parsing eff.charges"
58+
for s in `seq 1 $EIG_NIONS` ; do
59+
ss=`printf "%03u" $s`
60+
awk '{print $2,$3,$4}' borncs.$ss | tail -3 > bornch.$ss
61+
done
62+
mkdir temps02 ; mv eigens.* borncs.* temps02/
63+
printf " ..done\n"
64+
65+
# parse matrices, multiply them and collect squares (giving intensities)
66+
printf "..multiplying matrices, summing "
67+
for nu in `seq 1 $EIG_NVIBS` ; do
68+
nuu=`printf "%03u" $nu`
69+
int=0.0
70+
for alpha in 1 2 3 ; do # summing over alpha coordinates
71+
sumpol=0.0
72+
for s in `seq 1 $EIG_NIONS` ; do # summing over atoms
73+
ss=`printf "%03u" $s`
74+
awk -v a="$alpha" '(NR==a){print}' bornch.$ss > z.ion.$ss.alpha.$alpha
75+
# summing over beta coordinates and multiplying Z(s,alpha)*e(s) done by the following awk script
76+
paste z.ion.$ss.alpha.$alpha e.vib.$nuu.ion.$ss | \
77+
awk '{pol=$1*$4+$2*$5+$3*$6; print $0," ",pol}' > matr-vib-${nuu}-alpha-${alpha}-ion-${ss}
78+
done
79+
sumpol=`cat matr-vib-${nuu}-alpha-${alpha}-ion-* | awk '{sum+=$7} END {print sum}'`
80+
int=`echo "$int+($sumpol)^2" | sed 's/[eE]/*10^/g' | bc -l`
81+
done
82+
freq=`awk '(NR==1){print $8}' temps02/eigens.vib.$nuu`
83+
echo "$nuu $freq $int">> exact.res.txt
84+
printf "."
85+
done
86+
printf " ..done\n"
87+
88+
# format results, normalize intensities
89+
printf "..normalizing intensities"
90+
max=`awk '(NR==1){max=$3} $3>=max {max=$3} END {print max}' exact.res.txt`
91+
awk -v max="$max" '{printf "%03u %6.1f %5.3f\n",$1,$2,$3/max}' exact.res.txt > results.txt
92+
printf " ..done\n"
93+
94+
# clean up, display results
95+
printf "..finalizing:\n"
96+
mkdir temps03; mv bornch.* e.vib.*.allions temps03/
97+
mkdir temps04; mv z.ion* e.vib.*.ion.* temps04/
98+
mkdir temps05; mv matr-* temps05/
99+
mkdir results; mv *res*txt results/
100+
let NMATRIX=$EIG_NVIBS**2
101+
printf "%5u atoms found\n%5u vibrations found\n%5u matrices evaluated" \
102+
$EIG_NIONS $EIG_NVIBS $NMATRIX > results/statistics.txt
103+
# fast switch to clean up all temporary files
104+
rm -r temps*
105+
cat results/results.txt

‎VASP.py

+13-23
Original file line numberDiff line numberDiff line change
@@ -2,36 +2,26 @@
22
import subprocess
33
# 注意：所有函数读取的坐标统一为笛卡尔坐标
44

5-
def read_total_atoms():
6-
pipe = subprocess.Popen("sed -n '7p' POSCAR", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
7-
(content, error) = (pipe.stdout.read().strip(), pipe.stderr.read())
8-
if content != "":
9-
space = re.compile(r'\s+')
10-
content = space.split(content)
11-
number = 0
12-
for num in content:
13-
number += int(num)
14-
return number
155

16-
if error != "":
17-
print('')
18-
print('POSCAR: No such file or incorrect')
19-
number = input("please input the total number of atoms:")
20-
return int(number)
6+
class CmdRrror(Exception):
7+
def __init__(self, errorinfo):
8+
super().__init__(self) # 初始化父类
9+
self.errorinfo = errorinfo
10+
11+
def __str__(self):
12+
return 'Command Execution Error: ' + self.errorinfo
2113

2214

23-
def grep_outcar(command):
15+
def execCmd(command):
2416
pipe = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
25-
(content, error) = (pipe.stdout.readlines(), pipe.stderr.read())
17+
(content, error) = (pipe.stdout.readlines(), pipe.stderr.read().decode())
2618
if content:
19+
for i in range(len(content)):
20+
content[i] = content[i].decode()
2721
return content
2822

2923
if error != "":
30-
print('')
31-
print('POSCAR: No such file')
32-
print('')
33-
exit(0)
34-
return []
24+
raise CmdRrror(error)
3525

3626

3727
# this function can only read file in VASP 5.0 or later
@@ -183,7 +173,7 @@ def readGjf(file_name):
183173
return elements[1:], num_atoms, coordinates
184174

185175

186-
def writeGif(file_name, elements, num_atoms, coordinates):
176+
def writeGjf(file_name, elements, num_atoms, coordinates):
187177
"""
188178
:param file_name: str
189179
:param elements: [element1, element2, ...], str

‎cel2vas.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -18,4 +18,4 @@
1818
writeVasp(file_name.replace('.cell', '.vasp'), 1.0, basis, elements,
1919
num_atoms, '', coordinate_type, coordinates, [])
2020

21-
print(' ----------------- Done -----------------')
21+
print(' ----------------- Done -----------------\n')

‎chgflag.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@
1818
flag = sys.argv[-2]
1919
if flag != 'T' and flag != 'F':
2020
print('')
21-
print("error: unidentified flag %s" % flag)
21+
print("Error: unidentified flag '%s'" % flag)
2222
print('')
2323
exit(1)
2424

‎drawcluster.py

+92
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
#!/bin/env python3
2+
3+
import sys
4+
import os
5+
from VASP import readCell
6+
import re
7+
import numpy as np
8+
9+
"""
10+
crystal_conf.txt 格式:
11+
x, y, z //在x,y,z三个方向上晶胞的数量
12+
h, k, l, distance //晶面指标,距离原点的距离
13+
...
14+
h, k, l, distance //晶面指标,距离原点的距离
15+
"""
16+
17+
if len(sys.argv) == 1:
18+
print('')
19+
print('Usage: %s cell_file' % sys.argv[0].split('/')[-1])
20+
print('')
21+
exit(1)
22+
23+
if not os.path.isfile('crystal_conf.txt'):
24+
print('')
25+
print('Error: crystal_conf.txt missing!')
26+
print('')
27+
exit(1)
28+
29+
30+
def product(coord, hkl):
31+
# 截距为0时该项置0
32+
if hkl == 0:
33+
return np.zeros(coord.shape[0])
34+
else:
35+
return coord / hkl
36+
37+
epsilion = 0.000001
38+
space = re.compile(r'\s+')
39+
with open('crystal_conf.txt', 'r') as conf_file:
40+
content = conf_file.readlines()
41+
x, y, z = map(int, content[0].strip().split(','))
42+
43+
planes = []
44+
for i in range(1, len(content)):
45+
line = content[i].strip()
46+
if line == '':
47+
break
48+
49+
if line.startswith('#'):
50+
continue
51+
line = line.split(',')
52+
planes.append([int(line[0]), int(line[1]), int(line[2]), float(line[3])])
53+
54+
basis, elements, num_atoms, coordinate_type, coordinates = readCell(sys.argv[1])
55+
basis = np.array(basis)
56+
coordinates = np.array(coordinates)
57+
58+
# 构建超胞
59+
coordinates_ret = []
60+
for ix in range(-x, x + 1):
61+
coordinates_ret.append(coordinates + basis[0] * ix)
62+
coordinates = np.concatenate(coordinates_ret)
63+
64+
coordinates_ret = []
65+
for iy in range(-y, y + 1):
66+
coordinates_ret.append(coordinates + basis[1] * iy)
67+
coordinates = np.concatenate(coordinates_ret)
68+
69+
coordinates_ret = []
70+
for iz in range(-z, z + 1):
71+
coordinates_ret.append(coordinates + basis[2] * iz)
72+
coordinates = np.concatenate(coordinates_ret)
73+
74+
print('')
75+
for i, plane in enumerate(planes):
76+
# 知道截距时候的平面公式: x/h + y/k + z/l = 1
77+
print('processing plane: (%2d, %2d, %2d); distance = %5.4f' % (plane[0], plane[1], plane[2], plane[3]))
78+
ret = np.zeros(coordinates.shape[0])
79+
ret += product(coordinates[:, 0], plane[0] * plane[3])
80+
ret += product(coordinates[:, 1], plane[1] * plane[3])
81+
ret += product(coordinates[:, 2], plane[2] * plane[3])
82+
83+
coordinates = coordinates[ret < 1.0]
84+
print('')
85+
86+
with open('ret.xyz', 'w') as output_file:
87+
output_file.write("%d\n" % coordinates.shape[0])
88+
output_file.write('create from python\n')
89+
for coord in coordinates:
90+
output_file.write("%2s %16.10f %16.10f %16.10f\n" % (elements[0], coord[0], coord[1], coord[2]))
91+
output_file.write('\n')
92+

‎excongjf.py

+8-7
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919

2020
print('')
2121
print('--------------------')
22+
print('')
2223
print("Default file: %s" % default_file)
2324
dir_list = []
2425
for directory in sys.argv[1:]:
@@ -30,13 +31,13 @@
3031
print('')
3132
print('Processing...')
3233
for directory in dir_list:
33-
pipe = subprocess.Popen("vas2gjf.py %s" % (directory + '/' + default_file),
34-
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
35-
error = pipe.stderr.read()
36-
if error != '':
37-
print("Error in path(%s): %s" % (directory, error))
38-
continue
39-
os.system("mv %s %s" % ('%s/%s.gjf' % (directory, default_file), '%s.gjf' % directory))
34+
pipe = subprocess.Popen("vas2gjf.py %s" % (directory + '/' + default_file),
35+
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
36+
error = pipe.stderr.read().decode()
37+
if error != '':
38+
print("Error in path(%s): %s" % (directory, error))
39+
continue
40+
os.system("mv %s/%s.gjf %s" % (directory, default_file, '%s.gjf' % directory))
4041
print('')
4142
print('----- Done -----')
4243
print('')

0 commit comments

Comments
 (0)