forked from WMD-group/MacroDensity
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathInputDifference.py
executable file
·90 lines (82 loc) · 4.13 KB
/
InputDifference.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
#! /usr/bin/env python
import NewPotentialModule as pot
import math
import numpy as np
import matplotlib.pyplot as plt
#------------------------------------------------------------------
# READING
# Get the two potentials and change them to a planar average.
# This section should not be altered
#------------------------------------------------------------------
# SLAB
vasp_pot, NGX, NGY, NGZ, Lattice = pot.read_vasp_density('CHGCAR.Slab')
mag_a,mag_b,mag_c,vec_a,vec_b,vec_c = pot.matrix_2_abc(Lattice)
resolution_x = mag_a/NGX
resolution_y = mag_b/NGY
resolution_z = mag_c/NGZ
Volume = pot.get_volume(vec_a,vec_b,vec_c)
grid_pot_slab, electrons_slab = pot.density_2_grid(vasp_pot,NGX,NGY,NGZ,True,Volume)
# Save the lattce vectors for use later
Vector_A = [vec_a,vec_b,vec_c]
#----------------------------------------------------------------------------------
# CONVERT TO PLANAR DENSITIES
#----------------------------------------------------------------------------------
planar_slab = pot.planar_average(grid_pot_slab,NGX,NGY,NGZ)
# BULK
vasp_pot, NGX, NGY, NGZ, Lattice = pot.read_vasp_density('CHGCAR.Bulk')
mag_a,mag_b,mag_c,vec_a,vec_b,vec_c = pot.matrix_2_abc(Lattice)
resolution_x = mag_a/NGX
resolution_y = mag_b/NGY
resolution_z = mag_c/NGZ
# Save the lattce vectors for use later
Vector_B = [vec_a,vec_b,vec_c]
Volume = pot.get_volume(vec_a,vec_b,vec_c)
#----------------------------------------------------------------------------------
# CONVERT TO PLANAR DENSITIES
#----------------------------------------------------------------------------------
grid_pot_bulk, electrons_bulk = pot.density_2_grid(vasp_pot,NGX,NGY,NGZ,True,Volume)
planar_bulk = pot.planar_average(grid_pot_bulk,NGX,NGY,NGZ)
#----------------------------------------------------------------------------------
# FINISHED READING
#----------------------------------------------------------------------------------
# GET RATIO OF NUMBERS OF ELECTRONS
#----------------------------------------------------------------------------------
elect_ratio = int(electrons_slab/electrons_bulk)
#----------------------------------------------------------------------------------
# SPLINE THE TWO GENERATING A DISTANCE ON ABSCISSA
#----------------------------------------------------------------------------------
slab, bulk = pot.matched_spline_generate(planar_slab,planar_bulk,Vector_A[2],Vector_B[1])
#----------------------------------------------------------------------------------
# EXTEND THE BULK POTENTIAL TO MATCH THE ELECTRON NUMBERS
#----------------------------------------------------------------------------------
bulk = pot.extend_potential(bulk,elect_ratio,Vector_B[1])
#----------------------------------------------------------------------------------
# MATCH THE RESOLUTIONS OF THE TWO
#----------------------------------------------------------------------------------
slab, bulk = pot.match_resolution(slab, bulk)
plt.plot(bulk[:,1])
plt.show()
#----------------------------------------------------------------------------------
# TRANSLATE THE BULK POTENTIAL TO GET OVERLAP
#----------------------------------------------------------------------------------
bulk_trans = pot.translate_grid(bulk, 3.13,True,np.dot(Vector_B[1],elect_ratio),0.42)
bulk_trans = pot.translate_grid(bulk_trans, 6.57,False,np.dot(Vector_B[1],elect_ratio),0.42)
slab_trans = pot.translate_grid(slab, 6.5653,True,Vector_A[2])
#potential_difference = pot.diff_potentials(pot_slab_orig,bulk_extd,10,40,tol=0.04)
#plt.plot(bulk[:,0],bulk[:,1])
## PLOT THE POTENTIALS, THIS IS A CHECK THAT THEY HAVE ALIGNED CORRECTLY
#plt.plot(bulk_trans[:,0],bulk_trans[:,1],)
#plt.plot(slab_trans[:,0],slab_trans[:,1],)
#plt.show()
##------------------------------------------------------------------
## SET THE CHARGE DENSITY TO ZERO OUTSIDE THE BULK
bulk_vacuum = pot.bulk_vac(bulk_trans, slab_trans)
plt.plot(bulk_vacuum[:,0],bulk_vacuum[:,1])
plt.plot(slab_trans[:,0],slab_trans[:,1],)
plt.show()
# GET THE DIFFERENCES (within a numerical tolerence)
difference = pot.subs_potentials(slab_trans,bulk_vacuum,tol=0.01)
difference = pot.spline_generate(difference)
plt.plot(difference[:,0],difference[:,1])
#plt.plot(difference[:,0],difference[:,1])
plt.show()