#!/usr/bin/env python ######################################### import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt # import numpy as np # from distribution import gauss_bin # from matplotlib import rc # ######################################### rc('mathtext', fontset ='stix') #angle = 0 #-------------------------------------------------------------- # Preparing arrays #-------------------------------------------------------------- surf_atoms = [[ 2.01008073e+00, 1.42134188e+00, 1.11529848e+01], [ 2.01008073e+00, 4.26402539e+00, 1.11529848e+01], [ 2.01008073e+00, 7.10670890e+00, 1.11529848e+01], [ -2.17850106e-02, 1.17300492e-07, 9.99517836e+00], [ -2.17850106e-02, 2.84268362e+00, 9.99517836e+00], [ -2.17850106e-02, 5.68536713e+00, 9.99517836e+00], [ 4.04194649e+00, 1.15560770e-07, 9.99517818e+00], [ 4.04194649e+00, 2.84268362e+00, 9.99517818e+00], [ 4.04194649e+00, 5.68536713e+00, 9.99517818e+00], # [ 8.01853813e+00, 1.17300492e-07, 9.99517836e+00], [ 8.01853813e+00, 2.84268362e+00, 9.99517836e+00], [ 8.01853813e+00, 5.68536713e+00, 9.99517836e+00], [ -2.17850106e-02, 8.52805064e+00, 9.99517836e+00], [ 4.04194649e+00, 8.52805063e+00, 9.99517818e+00], [ 8.01853813e+00, 8.52805064e+00, 9.99517836e+00], # [ 6.03024227e+00, 1.42134184e+00, 8.74980064e+00], [ 6.03024227e+00, 4.26402535e+00, 8.74980064e+00], [ 6.03024227e+00, 7.10670885e+00, 8.74980064e+00], [ 2.01008070e+00, 1.42134182e+00, 8.40328480e+00], [ 2.01008070e+00, 4.26402533e+00, 8.40328480e+00], [ 2.01008070e+00, 7.10670883e+00, 8.40328480e+00], [ -7.06894583e-02, 4.51048592e-08, 7.13999074e+00], [ -7.06894583e-02, 2.84268355e+00, 7.13999074e+00], [ -7.06894583e-02, 5.68536706e+00, 7.13999074e+00], [ 4.09085093e+00, 4.45505359e-08, 7.13999066e+00], [ 4.09085093e+00, 2.84268355e+00, 7.13999066e+00], [ 4.09085093e+00, 5.68536706e+00, 7.13999066e+00], # DOESN'T INCLUDE THERMAL EXPANSION COEFFICIENT # [ 7.96963368e+00, 4.51048592e-08, 7.13999074e+00]] # xcell = 8.0830296493066 ycell = 8.5733476163356 zcell = 24.2122243358601 init_COM = [] close_COM = [] rxn_COM_D = [] rxn_COM_H = [] trap_COM = [] for i in range (len(surf_atoms)): for j in range (3): surf_atoms[i][j] = surf_atoms[i][j] * 1.005311542 # Thermal expansion coefficient #--------------------------------------------------------- # Reading data #--------------------------------------------------------- fdata = open('initial_COM.dat','r') # initial x, initial y, initial z data_init = fdata.readlines() fdata.close() for i in range (len(data_init)): COM=([float(x) for x in data_init[i].split()]) init_COM.append(COM) fdata = open('react_COM_D.dat','r') # initial x, initial y, initial z, react x, react y, react z data_react_D = fdata.readlines() fdata.close() for i in range (len(data_react_D)): COM_D=([float(x) for x in data_react_D[i].split()]) if COM_D[3] < 0: COM_D[3] += xcell if COM_D[3] > xcell: COM_D[3] -= xcell if COM_D[4] < 0: COM_D[4] += ycell if COM_D[4] > ycell: COM_D[4] -= ycell if COM_D[0] < 0: COM_D[0] += xcell if COM_D[0] > xcell: COM_D[0] -= xcell if COM_D[1] < 0: COM_D[1] += ycell if COM_D[1] > ycell: COM_D[1] -= ycell rxn_COM_D.append(COM_D) fdata = open('react_COM_H.dat','r') # initial x, initial y, initial z, react x, react y, react z data_react_H = fdata.readlines() fdata.close() for i in range (len(data_react_H)): COM_H=([float(x) for x in data_react_H[i].split()]) if COM_H[3] < 0: COM_H[3] += xcell if COM_H[3] > xcell: COM_H[3] -= xcell if COM_H[4] < 0: COM_H[4] += ycell if COM_H[4] > ycell: COM_H[4] -= ycell if COM_H[0] < 0: COM_H[0] += xcell if COM_H[0] > xcell: COM_H[0] -= xcell if COM_H[1] < 0: COM_H[1] += ycell if COM_H[1] > ycell: COM_H[1] -= ycell rxn_COM_H.append(COM_H) fdata = open('trap_COM.dat','r') # initial x, initial y, initial z, react x, react y, react z data_trap = fdata.readlines() fdata.close() for i in range (len(data_trap)): COM_trp=([float(x) for x in data_trap[i].split()]) if COM_trp[3] < 0: COM_trp[3] += xcell if COM_trp[3] > xcell: COM_trp[3] -= xcell if COM_trp[4] < 0: COM_trp[4] += ycell if COM_trp[4] > ycell: COM_trp[4] -= ycell if COM_trp[0] < 0: COM_trp[0] += xcell if COM_trp[0] > xcell: COM_trp[0] -= xcell if COM_trp[1] < 0: COM_trp[1] += ycell if COM_trp[1] > ycell: COM_trp[1] -= ycell trap_COM.append(COM_trp) #------------------------------------------------------------------------ # Playing with arrays #------------------------------------------------------------------------ surf_atoms=np.array(surf_atoms) surf_atoms=surf_atoms.transpose() init_COM=np.array(init_COM) init_COM=init_COM.transpose() rxn_COM_D=np.array(rxn_COM_D) rxn_COM_D=rxn_COM_D.transpose() rxn_COM_H=np.array(rxn_COM_H) rxn_COM_H=rxn_COM_H.transpose() trap_COM=np.array(trap_COM) trap_COM=trap_COM.transpose() #------------------------------------------------------------------------ # Making plots - points #------------------------------------------------------------------------ # xy plot plt.figure(figsize=(xcell+2,ycell+2)) plt.plot(surf_atoms[0,0:3], surf_atoms[1,0:3], color='lightgrey', mew = 2.5, marker='o', markersize = 100, linestyle='None') plt.plot(surf_atoms[0,3:15], surf_atoms[1,3:15], color='lightgrey', mew = 1.5, marker='o', markersize = 100, linestyle='None') plt.plot(surf_atoms[0,15:18], surf_atoms[1,15:18], color='lightgrey', mew = 0.5, marker='o', markersize = 100, linestyle='None') #plt.plot(init_xy[14],init_xy[15],'ro') plt.plot(init_COM[0], init_COM[1], 'wo', mec = 'k', markersize = 10) plt.plot(rxn_COM_D[0], rxn_COM_D[1], 'wo', mec = 'b', markersize = 10, mew = 1) plt.plot(rxn_COM_H[0], rxn_COM_H[1], 'wo', mec = 'r', markersize = 10, mew = 1) plt.plot(trap_COM[0], trap_COM[1], 'wo', mec = 'g', markersize = 10, mew = 1) plt.plot(rxn_COM_D[0], rxn_COM_D[1], 'x', mec = 'b', markersize = 8, mew = 1) plt.plot(rxn_COM_H[0], rxn_COM_H[1], 'x', mec = 'r', markersize = 8, mew = 1) plt.plot(trap_COM[0], trap_COM[1], 'x', mec = 'g', markersize = 8, mew = 1) plt.plot(rxn_COM_D[3], rxn_COM_D[4], 'bo', markersize = 10) plt.plot(rxn_COM_H[3], rxn_COM_H[4], 'ro',markersize = 10) plt.plot(trap_COM[3], trap_COM[4], 'go', markersize = 10) plt.xlim(-1,xcell+1) plt.ylim(-1,ycell+1) plt.xticks(()) plt.yticks(()) plt.plot(surf_atoms[0,0:3], surf_atoms[1,0:3], color='None', mew = 2.5, marker='o', markersize = 100, linestyle='None') plt.plot(surf_atoms[0,3:15], surf_atoms[1,3:15], color='None', mew = 1.5, marker='o', markersize = 100, linestyle='None') plt.plot(surf_atoms[0,15:18], surf_atoms[1,15:18], color='None', mew = 0.5, marker='o', markersize = 100, linestyle='None') #plt.savefig('xy_points_python.eps',bbox_inches='tight') plt.savefig('Fig6.png', dpi = 600) #plt.show()