import pyiron, ase.io
import numpy as np
import matplotlib.pyplot as plt
import os.path
from ase.io import lammpsdata
%config InlineBackend.figure_format = 'retina'
pr = pyiron.Project('ADIS')
bulk = pr.create_ase_bulk('Cu', cubic=True).repeat(10)
vacancy = bulk.copy()
vacancy_id = 17
vacancy[vacancy_id]
vac_position = vacancy[vacancy_id].position
del vacancy[vacancy_id]
vacancy.plot3d()
NGLWidget()
voro = vacancy.analyse_ovito_voronoi_volume()
plt.xlabel('Voronoi volume')
plt.hist(voro, log=True);
vac_neighbors = voro > 12
vac_neighbors.sum()
12
vacancy[vac_neighbors].plot3d()
NGLWidget()
vacancy_neighbors_structure = vacancy[vac_neighbors]
neighbors = vacancy_neighbors_structure.get_neighbors(num_neighbors=11)
print('Predicted position:',
np.round(vacancy_neighbors_structure.positions[0]+neighbors.vecs[0, :, :].sum(axis = 0)/12, 8))
print('Original position:',
np.round(vac_position, 8))
Predicted position: [ 0. 16.245 1.805] Original position: [ 0. 16.245 1.805]
cu_400_5 = pr.load('Cu_400_5').structure
cu_400_5[:] = 'Cu'
cu_400_5.plot3d()
NGLWidget()
voro = cu_400_5.analyse_ovito_voronoi_volume()
plt.xlabel('Voronoi volume')
plt.hist(voro, log=True, bins=50);
j = pr.create_job(pr.job_type.Lammps, "Cu_minimization", True)
j.structure = cu_400_5
j.potential = j.list_potentials()[0]
j.calc_minimize()
j.run()
The job Cu_minimization was saved and received the ID: 48
j.animate_structure()
NGLWidget(max_frame=1)
cu_400_5_q = j.get_structure()
voro = cu_400_5_q.analyse_ovito_voronoi_volume()
plt.hist(voro, bins=100, log=True);
(voro > 12).sum()
58
vacancy_neighbors_structure = cu_400_5_q[voro>12]
vacancy_neighbors_structure.center_coordinates_in_unit_cell();
vacancy_neighbors_structure.plot3d()
NGLWidget()
neighbors = vacancy_neighbors_structure.get_neighbors(num_neighbors=11)
plt.xlabel('Distance')
plt.hist(neighbors.distances.max(axis=-1));
vacancy_neighbors_structure[neighbors.distances.max(axis=-1) < 4] = 'Ni'
vacancy_neighbors_structure.plot3d()
NGLWidget()
vac_positions = vacancy_neighbors_structure.positions + neighbors.vecs.sum(axis = 1)/12
vac_structure = pr.create_atoms(positions = vac_positions, symbols = ['O']*len(vac_positions), pbc=True,
cell = vacancy_neighbors_structure.cell)
vac_structure.append(vacancy_neighbors_structure)
vac_structure.plot3d()
NGLWidget()
import scipy.stats as st
P = vacancy.get_extended_positions(10)[0].T
kde = st.gaussian_kde(P, bw_method=.15)
x = np.linspace(0, 36.1, 50)
X, Y, Z = np.meshgrid(x, x, x)
positions = np.vstack([X.ravel(), Y.ravel(), Z.ravel()])
image = np.reshape(kde(positions).T, X.shape)
positions
array([[ 0. , 0. , 0. , ..., 36.1 , 36.1 , 36.1 ], [ 0. , 0. , 0. , ..., 36.1 , 36.1 , 36.1 ], [ 0. , 0.73673469, 1.47346939, ..., 34.62653061, 35.36326531, 36.1 ]])
W = np.where(Z < .75e-5)
plt.style.use('default')
# plt.rc('figure', figsize=(14,10))
np.unique(P[2][P[2]>0])
array([ 1.805, 3.61 , 5.415, 7.22 , 9.025, 10.83 , 12.635, 14.44 , 16.245, 18.05 , 19.855, 21.66 , 23.465, 25.27 , 27.075, 28.88 , 30.685, 32.49 , 34.295, 36.1 , 37.905, 39.71 , 41.515, 43.32 , 45.125])
# plane = P[2, 2000]
plane = np.unique(P[2][P[2]>0])[10]
z_index = np.argmin(np.absolute(x-plane))
# print(plane)
plt.imshow(np.rot90(np.rot90(image[:, :, z_index])), extent=(0, 36.1, 0, 36.1),)
plt.colorbar()
# plt.scatter(P[0, abs(P[2] - plane) < 0.5], P[1, abs(P[2] - plane) < 0.5], marker='.', color='r')
plt.scatter(P[0, abs(P[2] - plane) < 0.5], P[1, abs(P[2] - plane) < 0.5], marker='.', color='r')
plt.xlim(0, 36.1)
plt.ylim(0, 36.1);
plt.xlabel('Gaussian density')
# plt.axvline(5.7055e-6, color='r')
plt.hist(image.flat, log=True, bins=100);