-
Notifications
You must be signed in to change notification settings - Fork 1
/
MANULS_plot_1D.py
115 lines (69 loc) · 2.83 KB
/
MANULS_plot_1D.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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
font = {'family' : 'normal',
'weight' : 'normal',
'size' : 14}
matplotlib.rc('font', **font)
######################################################################
##################### USER DEFINED SECTION ##########################
######################################################################
# This is the only section that needs to be modified, the rest is up to the code
# Loading of the files containing the data, please see the documentation
data=np.loadtxt('/path/to/file_grid.txt')
polar=np.loadtxt('/path/to/file_polarizability.txt')
# "pointsx" and "pointsy" = number of points scanned for each variable
pointsx=int(37)
#the optimal electric field
#e=np.array([field_x, field_y, field_z])
e=np.array([-0.00395539, 0.00419754, 0.00640764])
######################################################################
##################### END OF USER DEFINED SECTION ###################
######################################################################
########################################################################
########################################################################
########################################################################
energy=data[:,1]
dipx=data[:,2]
dipy=data[:,3]
dipz=data[:,4]
energy_column=data[:,1]
a00=np.array([])
a01=np.array([])
a02=np.array([])
a10=np.array([])
a11=np.array([])
a12=np.array([])
a20=np.array([])
a21=np.array([])
a22=np.array([])
for i in range(int(np.shape(polar)[0]/3)):
a00=np.append(a00,polar[0+i*3,0])
a01=np.append(a01,polar[0+i*3,1])
a02=np.append(a02,polar[0+i*3,2])
a10=np.append(a10,polar[1+i*3,0])
a11=np.append(a11,polar[1+i*3,1])
a12=np.append(a12,polar[1+i*3,2])
a20=np.append(a20,polar[2+i*3,0])
a21=np.append(a21,polar[2+i*3,1])
a22=np.append(a22,polar[2+i*3,2])
fig=plt.figure()
x=np.linspace(data[0,0],data[-1,0],num=pointsx)
plt.plot(x,(energy-np.amin(energy))*627.503,marker='o')
plt.title('Unperturbed PES')
plt.xlabel(r'Variable Scanned')
plt.ylabel(r'Relative Energy (kcal/mol) ')
plt.show()
#####POLARIZABILITY
energy_perturbed=np.zeros(pointsx)
for g in range(int(pointsx)):
energy_perturbed[g]=energy[g]-(dipx[g]*e[0]+dipy[g]*e[1]+dipz[g]*e[2])-(1/2)*e[0]*(e[0]*a00[g]+e[1]*a01[g]+e[2]*a02[g])-(1/2)*e[1]*(e[0]*a10[g]+e[1]*a11[g]+e[2]*a12[g])-(1/2)*e[2]*(e[0]*a20[g]+e[1]*a21[g]+e[2]*a22[g])
fig=plt.figure()
x=np.linspace(data[0,0],data[-1,0],num=pointsx)
plt.plot(x,(energy_perturbed-np.amin(energy))*627.503,marker='o',color='tab:orange')
plt.title('Perturbed PES')
plt.xlabel(r'Variable Scanned')
plt.ylabel(r'Relative Energy (kcal/mol)')
plt.show()