forked from WMD-group/MacroDensity
-
Notifications
You must be signed in to change notification settings - Fork 0
/
FieldAtPoint.py
executable file
·130 lines (104 loc) · 4.8 KB
/
FieldAtPoint.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
124
125
126
127
128
129
#! /usr/bin/env python
# FieldAtPoint.py - try to calculate the electric field (grad of potential) at arbitrary point
# Forked form PlaneField.py - JMF 2016-01-25
import macrodensity as md
import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors,cm #colour maps; so I can specify cube helix
import sys #for argv
## Input section (define the plane with 3 points, fractional coordinates)
a_point = [0, 0, 0]
b_point = [1, 0, 1]
c_point = [0, 1, 0]
#LOCPOT.CsPbI3_cubic LOCPOT.CsPbI3_distorted LOCPOT.MAPI_pseudocubic
#input_file = 'LOCPOT.CsPbI3_distorted'
input_file = sys.argv[1]
print("Input file ",input_file)
#------------------------------------------------------------------
# Get the potential
# This section should not be altered
#------------------------------------------------------------------
vasp_pot, NGX, NGY, NGZ, Lattice = md.read_vasp_density(input_file)
vector_a,vector_b,vector_c,av,bv,cv = md.matrix_2_abc(Lattice)
resolution_x = vector_a/NGX
resolution_y = vector_b/NGY
resolution_z = vector_c/NGZ
grid_pot, electrons = md.density_2_grid(vasp_pot,NGX,NGY,NGZ)
## Get the gradiens (Field), if required.
## Comment out if not required, due to compuational expense.
print("Calculating gradients (Electic field, E=-Grad.V )...")
grad_x,grad_y,grad_z = np.gradient(grid_pot[:,:,:],resolution_x,resolution_y,resolution_z)
#------------------------------------------------------------------
##------------------------------------------------------------------
## Get the equation for the plane
## This is the section for plotting on a user defined plane;
## uncomment commands if this is the option that you want.
##------------------------------------------------------------------
## Convert the fractional points to grid points on the density surface
a = md.numbers_2_grid(a_point,NGX,NGY,NGZ)
b = md.numbers_2_grid(b_point,NGX,NGY,NGZ)
c = md.numbers_2_grid(c_point,NGX,NGY,NGZ)
plane_coeff = md.points_2_plane(a,b,c)
## Calculate magnitude of gradient.
# Should be able to use numpy.linalg.norm for this, but the Python array indices are causing me grief
X2 = np.multiply(grad_x,grad_x)
Y2 = np.multiply(grad_y,grad_y)
Z2 = np.multiply(grad_z,grad_z)
grad_mag = np.sqrt(np.add(X2,Y2,Z2))
# This was my, non working, attempt to use the built in function.
#grad_mag=np.linalg.norm( [grad_y,grad_y,grad_z], axis=3)
## This function in Macrodensity averages Efield ACROSS Z for Slab calculations
#xx,yy,grd = pot.create_plotting_mesh(NGX,NGY,NGZ,plane_coeff,grad_mag) #AVG over full volume
# Here we construct the same xx,yy,grd variables with a SLICE, forming a plane in XY at particular ZSLICE
xx, yy = np.mgrid[0:NGX,0:NGY]
ZSLICE= NGZ/2 # Chosses where in the Z axis the XY slice is cut through
# Slice of magnitude of electric field, for contour plotting
grd=grad_mag[:,:,ZSLICE]
# Slices of x and y components for arrow plotting
grad_x_slice=grad_x[:,:,ZSLICE]
grad_y_slice=grad_y[:,:,ZSLICE]
# OK, that's all our data
# This code tiles the data to (2,2) to re-expand unit cell to a 2x2 supercell in XY
xx,yy=np.mgrid[0:2*NGX,0:2*NGY]
grd=np.tile(grd, (2,2))
grad_x_slice=np.tile(grad_x_slice, (2,2))
grad_y_slice=np.tile(grad_y_slice, (2,2))
# End of tiling code
## Contours of the above sliced data
plt.contour(xx,yy,grd,6,cmap=cm.cubehelix)
# Also generate a set of Efield arrows ('quiver') for this data.
# Specifying the drawing parameters is quite frustrating - they are very brittle + poorly documented.
plt.quiver(xx,yy, grad_x_slice, grad_y_slice,
color='grey',
units='dots', width=1, headwidth=3, headlength=4
) #,
# units='xy', scale=10., zorder=3, color='blue',
# width=0.007, headwidth=3., headlength=4.)
plt.axis('equal') #force square aspect ratio; this assuming X and Y are equal.
plt.show()
##------------------------------------------------------------------
##------------------------------------------------------------------
# CsPbI3 - distorted
PB_X=0.469972*NGX
PB_Y=0.530081*NGY
PB_Z=0.468559*NGZ
# CsPbI3 - perfect cubic
#PB_X=0.5*NGX
#PB_Y=0.5*NGY
#PB_Z=0.5*NGZ
# MAPBI3 - pseudo cubic distorted
#PB_X=0.476171*NGX
#PB_Y=0.500031*NGY
#PB_Z=0.475647*NGZ
# Read out massive grad table, in {x,y,z} components
print(grad_x[PB_X][PB_Y][PB_Z],grad_y[PB_X][PB_Y][PB_Z],grad_z[PB_X][PB_Y][PB_Z])
# Norm of electric field at this point
print(np.linalg.norm([ grad_x[PB_X][PB_Y][PB_Z],grad_y[PB_X][PB_Y][PB_Z],grad_z[PB_X][PB_Y][PB_Z] ]))
# OK, let's try this with a Spectral method (FFT)
# JMF - Not currently working; unsure of data formats, need worked example
from scipy import fftpack
V_FFT=fftpack.fftn(grid_pot[:,:,:])
V_deriv=fftpack.diff(grid_pot[:,:,:]) #V_FFT,order=1)
# Standard catch all to drop into ipython at end of script for variable inspection etc.
from IPython import embed; embed() # End on an interactive ipython console to inspect variables etc.