-
Notifications
You must be signed in to change notification settings - Fork 0
/
aims.v_eff.2.py
executable file
·127 lines (96 loc) · 3.29 KB
/
aims.v_eff.2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
#!/usr/bin/python
######################################
#
# Script which calculate average
# efective potential basing on the
# v_eff.dat file
#
######################################
import sys
import numpy as np
from ase.io.cube import write_cube
from ase.io.aims import read_aims
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
num = len(sys.argv)
f = open(sys.argv[num-1], "r")
lines = f.read().splitlines()
points_tmp = []
values_tmp = []
for l in lines:
r = l.split()
row = [float(r[0]), float(r[1]), float(r[2]), float(r[3])]
points_tmp.append(row[0:3])
values_tmp.append(row[3])
points3D = np.array(points_tmp)
values3D = np.array(values_tmp)
max_x = max(points3D[:,0])
max_y = max(points3D[:,1])
max_z = max(points3D[:,2])
min_x = min(points3D[:,0])
min_y = min(points3D[:,1])
min_z = min(points3D[:,2])
dd=int(sys.argv[num-2])
density=complex(0.0,dd)
npoints = density.imag
z0 = min_z + (max_z-min_z)*0.48
dz = ((max_z-min_z)/npoints)
#grid_x, grid_y = np.mgrid[min_x:max_x:2*density, min_y:max_y:density]
grid_x, grid_y = np.mgrid[min_x:max_x:density, min_y:max_y:2*density]
interpolation = np.empty([npoints*2,npoints])
#print "len(interpolation)", len(interpolation)
#print interpolation
#for z0 in np.arange(min_z,max_z,dz):
#print z0, min_z, max_z
#print np.arange(min_z,max_z,dz)
for z0 in np.arange(min_z,max_z,dz):
condition = (points3D[:,2] > z0) & (points3D[:,2] < z0+dz)
size=len(np.extract(condition, condition))
if(size>1):
# print z0, z0+dz
# print condition.any()
# print (points3D[:,2] > z0)
# print (points3D[:,2] < z0+dz)
# print condition
# print
# print "KK"
p2D = points3D[condition]
points2D = p2D[:,0:2]
values2D = values3D[condition]
# print "p" ,len(points2D), points2D
# print "v" ,len(values2D), values2D
# print "gx" ,len(grid_x), grid_x
# print "gy" ,len(grid_y), grid_y
#grid_x, grid_y , grid_z = np.mgrid[min_x:max_x:2*density, min_y:max_y:density, min_z:max_z:density]
#interpolation = griddata(points3D, values3D, (grid_x, grid_y , grid_z), method='nearest')
#interpolation = interpolation[:,:,npoints*0.47]
# print points2D
# print values2D
tmp = griddata(points2D, values2D, (grid_x, grid_y), method='nearest')
# print "t" ,len(tmp), tmp
plt.clf()
plt.grid()
# plt.axis([0.4*npoints,0.6*npoints,npoints,0.5*npoints])
plt.axis([0.2*npoints,npoints,0.2*npoints,0.8*npoints])
plt.imshow(tmp, interpolation='nearest')
plt.colorbar()
plt.savefig('v_eff_'+str((z0-min_z)/dz)+'j.png')
# plt.show()
# print tmp
# np.append(interpolation, tmp)
#print interpolation
#print len(interpolation)
#plt.imshow(interpolation, interpolation='nearest')
#plt.grid()
#plt.axis([0.4*npoints,0.6*npoints,0.5*npoints,npoints])
##plt.savefig('v_eff_'+str(npoints)+'j.png')
#plt.show()
#atoms = read_aims("geometry.in")
#cell = np.array( [[max_x-min_x, 0. ,0. ],\
# [0.0 , max_y-min_y,0. ],\
# [0. , 0. ,max_z-min_z]] )
#print cell
#atoms.set_cell(cell*0.5)
#shift = np.array([min_x,min_y,min_z])
#atoms.translate(-0.5*shift)
#write_cube("v_eff.cube",atoms,interpolation)