Skip to content

Commit 37db4d9

Browse files
qzhu2017qzhu2017
committed
remove get_rod/layer function
1 parent f046f00 commit 37db4d9

File tree

2 files changed

+111
-244
lines changed

2 files changed

+111
-244
lines changed

‎pyxtal/descriptor.py‎

Lines changed: 98 additions & 205 deletions
Original file line numberDiff line numberDiff line change
@@ -90,214 +90,107 @@ def get_q4_and_q6(center, neighbors, weights=None):
9090
if __name__ == '__main__':
9191

9292
from ase.build import bulk
93+
from ase import Atoms
9394
from ase.neighborlist import NeighborList
9495

95-
data = [
96-
('Si', 'diamond'),
97-
('Fe', 'bcc'),
98-
('Cu', 'fcc'),
99-
('Mg', 'hcp'),
100-
]
101-
102-
for d in data:
103-
(element, packing) = d
104-
crystal = bulk(element, packing)
96+
def calculate(element, packing, a, covera, cutoff):
97+
if packing in ['bct']:
98+
crystal = bulk(element, packing, a=a, c=a*covera)
99+
elif packing in ['tetragonal']:
100+
crystal = Atoms(element, cell=(a, a, a*covera), pbc=True,
101+
scaled_positions=[(0, 0, 0)])
102+
elif packing in ['diamond']:
103+
crystal = Atoms(8*element, cell=(a, a, a*covera), pbc=True,
104+
scaled_positions=[(0, 0, 0), (0.25, 0.25, 0.25),
105+
(0, 0.5, 0.5), (0.25, 0.75, 0.75),
106+
(0.5, 0, 0.5), (0.75, 0.25, 0.75),
107+
(0.5, 0.5, 0), (0.75, 0.75, 0.25)])
108+
else:
109+
crystal = bulk(element, packing, a=a, covera=covera)
105110
center = crystal.get_positions()[0]
106111
cell = crystal.get_cell()
107-
for cutoff in [2.4, 2.8, 3.0, 3.6, 4.2, 5.0]:
108-
rs = [cutoff/2]*len(crystal)
109-
nl = NeighborList(rs, self_interaction=False, bothways=True, skin=0.0)
110-
nl.update(crystal)
111-
indices, offsets = nl.get_neighbors(0)
112-
if len(indices) > 0:
113-
neighbors = np.zeros([len(indices),3])
114-
for i in range(len(indices)):
115-
pos, offset = crystal.positions[indices[i]], offsets[i]
116-
neighbors[i, :] = pos + np.dot(offset, cell)
117-
q4, q6 = get_q4_and_q6(center, neighbors, 'auto')
118-
print(element, packing, cutoff, len(neighbors), q4, q6)
119-
120-
121-
#import numpy as np
122-
#
123-
#class descriptor():
124-
# """
125-
# A class to handle the crystal descriptor
126-
#
127-
# Args:
128-
# xtal:
129-
# cutoff:
130-
# factor:
131-
# size:
132-
# """
133-
#
134-
# def __init__(self, struc, cutoff=5.0, factor=1.5, size=201):
135-
# self.struc = struc
136-
# self.cutoff = cutoff
137-
# self.factor = factor
138-
# self.size = size
139-
#
140-
# #initialize x, y
141-
# self.x = np.linspace(0, cutoff, size)
142-
# self.y = np.zeros([len(self.struc.mol_sites), size])
143-
# self.res = cutoff/(size-1)
144-
# self.compute()
145-
# #self.smear()
146-
#
147-
# def compute(self):
148-
# """
149-
# compute the intensities
150-
#
151-
# Args:
152-
# struc: pyxtal object
153-
# """
154-
# for i, site in enumerate(self.struc.mol_sites):
155-
# coord0, specie0 = site._get_coords_and_species(absolute=True, first=True)
156-
# res = self.struc.get_neighboring_molecules(i, self.factor, self.cutoff)
157-
# (min_ds, neighs, comps) = res
158-
# for neigh in neighs:
159-
# ds = cdist(coord0, neigh)
160-
# ds = ds[ds<self.cutoff]
161-
# d2 = ds**6
162-
# inv_a2 = 1/d2
163-
# intensities = inv_a2/np.sum(inv_a2) #intensity
164-
# for d, intensity in zip(ds, intensities):
165-
# j = round(d/self.res)
166-
# self.y[i, j] += intensity
167-
#
168-
# def plot(self, filename=None, figsize=(10,8), smears=[0.1, 0.5, 1.0]):
169-
#
170-
# import matplotlib.pyplot as plt
171-
#
172-
# plt.figure(figsize=figsize)
173-
#
174-
# x_min, x_max = 0, self.cutoff + 0.5
175-
# plt.gca()
176-
#
177-
# y_smears = []
178-
# for s in smears:
179-
# y_smears.append(self.smear(s))
180-
#
181-
# for id in range(self.y.shape[0]):
182-
# label = 'Site_{:d} (Raw)'.format(id)
183-
# plt.plot(self.x, self.y[id,:], label=label)
184-
# for j, s in enumerate(smears):
185-
# label = 'Site_{:d} ({:4.1f})'.format(id, s)
186-
# plt.plot(self.x, y_smears[j][id,:], label=label)
187-
# plt.legend()
188-
# plt.xlim([x_min, x_max])
189-
# plt.xlabel('R($\AA$)')
190-
# plt.ylabel('Intensity')
191-
#
192-
# if filename is None:
193-
# plt.show()
194-
# else:
195-
# plt.savefig(filename)
196-
# plt.close()
197-
#
198-
# def smear(self, sigma=0.8):
199-
# """
200-
# Apply Gaussian smearing to spectrum y value.
201-
# Args:
202-
# sigma: Std dev for Gaussian smear function
203-
# """
204-
# from scipy.ndimage.filters import gaussian_filter1d
205-
# y = np.zeros([len(self.struc.mol_sites), self.size])
206-
# for id in range(self.y.shape[0]):
207-
# y[id] = gaussian_filter1d(self.y[id], sigma)
208-
# return y
209-
#
210-
#
211-
#if __name__ == "__main__":
212-
#
213-
# from pyxtal import pyxtal
214-
# from pyxtal.XRD import Similarity
215-
# from pyxtal.representation import representation
216-
# from scipy.optimize import minimize
217-
# import warnings
218-
# warnings.filterwarnings("ignore")
219-
#
220-
#
221-
# #for code in ['ACSALA', 'ACSALA13']:
222-
# # c = pyxtal(molecular=True)
223-
# # c.from_CSD(code)
224-
# # des = descriptor(c)
225-
# # des.plot(code+'.png')
226-
# # print(np.sum(des.y))
227-
#
228-
# #c1 = pyxtal(molecular=True); c1.from_CSD('ACSALA'); des=descriptor(c1); f1=des.smear(1.0)
229-
# #c2 = pyxtal(molecular=True); c2.from_CSD('ACSALA13'); des=descriptor(c2); f2=des.smear(1.0)
230-
#
231-
# #smiles = ['CC(=O)Oc1ccccc1C(O)=O']
232-
# #rep1 = c1.get_1D_representation(); rep3 = representation(rep1.x, smiles); c3=rep3.to_pyxtal()
233-
# #rep2 = c2.get_1D_representation(); rep4 = representation(rep2.x, smiles); c4=rep4.to_pyxtal()
234-
#
235-
# #des=descriptor(c3); f3=des.smear(1.0)
236-
# #des=descriptor(c4); f4=des.smear(1.0)
237-
#
238-
# #f1 = [des.x, f1[0]]
239-
# #f2 = [des.x, f2[0]]
240-
# #f3 = [des.x, f3[0]]
241-
# #f4 = [des.x, f4[0]]
242-
#
243-
# #S12 = Similarity(f1, f2, l=0.1); print(S12.value)
244-
# #S13 = Similarity(f1, f3, l=0.1); print(S13.value)
245-
# #S23 = Similarity(f2, f3, l=0.1); print(S23.value)
246-
# #S24 = Similarity(f2, f4, l=0.1); print(S24.value)
247-
# #S34 = Similarity(f2, f4, l=0.1); print(S34.value)
248-
#
249-
# #
250-
# #c1.to_file('ACSALA.cif')
251-
# #c2.to_file('ACSALA13.cif')
252-
# #c3.to_file('rep-A.cif')
253-
# #c3.to_file('rep-B.cif')
254-
# ##print(S.value)
255-
# ##print(c)
256-
# ##print(c.mol_sites[0].molecule.smile)
257-
# ##rep = c.get_1D_representation()
258-
# #print(rep1.to_string())
259-
# #print(rep3.to_string())
260-
# #print(rep2.to_string())
261-
# #print(rep4.to_string())
262-
#
263-
# #print(c1)
264-
# #print(c3)
265-
# #print(c1.mol_sites[0].molecule.mol.cart_coords[:3,:])
266-
# #print(c3.mol_sites[0].molecule.mol.cart_coords[:3,:])
267-
# #d = np.zeros([len(f1[0]), 2])
268-
# #print(f1[0])
269-
# #d[:,0] = f1[0]
270-
# #d[:,1] = f1[1]
271-
# #np.savetxt('ref.txt', d)
272-
#
273-
# def fun(x, f0, smiles):
274-
# vec = [[14, 0, 11.45, 6.60, 11.39, 95.5], [0.22, 0.41, 0.03, x[0], x[1], x[2], -85.5, -2.4, -1.7, 1]]
275-
# rep = representation(vec, smiles)
276-
# c = rep.to_pyxtal()
277-
# des = descriptor(c)
278-
# f = des.smear(1.0)
279-
# f1 = [des.x, f[0]]
280-
# S = Similarity(f0, f1, l=0.1)
281-
# #print(x)
282-
# print(S.value)
283-
# return -S.value
284-
#
285-
# # 14 0 11.45 6.60 11.39 95.5 1 0.22 0.41 0.03 -42.7 26.4 32.0 -85.5 -2.4 -1.7 1
286-
# # 14 0 12.10 6.49 11.32 111.5 1 0.27 0.42 0.92 -135.9 26.1 -32.1 84.0 5.3 -1.5 1
287-
# smiles = ['CC(=O)Oc1ccccc1C(O)=O']
288-
# x0 = [-42.7, 26.4, 32.0]
289-
# f0 = np.loadtxt('ref.txt')
290-
# f0 = [f0[:,0], f0[:, 1]]
291-
# res = minimize(fun, x0, method='Powell', args=(f0, smiles), tol=1e-8, options={'maxiter':1000, 'disp': True})
292-
# #res = minimize(fun, x0, method='CG', args=(f0, smiles), tol=1e-6, options={'maxiter':1000, 'disp': True}) #opt w/o gradient
293-
# print(res.x)
294-
# x = res.x
295-
# vec = [[14, 0, 11.45, 6.60, 11.39, 95.5], [0.22, 0.41, 0.03, x[0], x[1], x[2], -85.5, -2.4, -1.7, 1]]
296-
# rep = representation(vec, smiles)
297-
# c = rep.to_pyxtal()
298-
# print(c)
299-
# c.to_file('ans.cif')
300-
#
301-
# #res = minimize(fun, x0, method='Powell', args=(f0), tol=1e-4, options={'maxiter':20}) #opt w/o gradient
112+
rs = [cutoff/2]*len(crystal)
113+
nl = NeighborList(rs, self_interaction=False, bothways=True, skin=0.0)
114+
nl.update(crystal)
115+
indices, offsets = nl.get_neighbors(0)
116+
if len(indices) > 0:
117+
neighbors = np.zeros([len(indices),3])
118+
for i in range(len(indices)):
119+
pos, offset = crystal.positions[indices[i]], offsets[i]
120+
neighbors[i, :] = pos + np.dot(offset, cell)
121+
q4, q6 = get_q4_and_q6(center, neighbors, 'auto')
122+
print(element, packing, cutoff, len(neighbors), q4, q6)
123+
124+
return q4, q6
125+
else:
126+
return None
127+
128+
_as = np.linspace(4.0, 8.0, 30) #[0, 4.0, 5.0, 6.0, 7.0]
129+
Npt = 80
130+
cut = _as[-1] + 1
131+
data = []
132+
for v in np.linspace(4, 16, Npt):
133+
data.append(('Xe', 'hcp', np.sqrt(v/3)))
134+
hcps = []
135+
for i, d in enumerate(data):
136+
(element, packing, covera) = d
137+
for j, a in enumerate(_as):
138+
q4, q6 = calculate(element, packing, a, covera, cut)
139+
hcps.append([q4, q6])
140+
141+
data = []
142+
for v in np.linspace(1, 8, Npt):
143+
data.append(('Xe', 'bct', np.sqrt(v/4)))
144+
bcts = []
145+
for i, d in enumerate(data):
146+
(element, packing, covera) = d
147+
for j, a in enumerate(_as):
148+
q4, q6 = calculate(element, packing, a, covera, cut)
149+
bcts.append([q4, q6])
150+
151+
data = []
152+
for v in np.linspace(1, 8, Npt):
153+
data.append(('Xe', 'tetragonal', np.sqrt(v/4)))
154+
ts = []
155+
for i, d in enumerate(data):
156+
(element, packing, covera) = d
157+
for j, a in enumerate(_as):
158+
q4, q6 = calculate(element, packing, a, covera, cut)
159+
ts.append([q4, q6])
160+
161+
data = []
162+
for v in np.linspace(1, 8, Npt):
163+
data.append(('Xe', 'diamond', np.sqrt(v/4)))
164+
165+
ds = []
166+
for i, d in enumerate(data):
167+
(element, packing, covera) = d
168+
for j, a in enumerate(_as+3):
169+
q4, q6 = calculate(element, packing, a, covera, cut)
170+
ds.append([q4, q6])
171+
172+
173+
hcps = np.array(hcps)
174+
bcts = np.array(bcts)
175+
ds = np.array(ds)
176+
ts = np.array(ts)
177+
import matplotlib.pyplot as plt
178+
plt.figure(figsize=[12, 12])
179+
plt.scatter(hcps[:,0], hcps[:,1], label='hcp')
180+
plt.scatter(bcts[:,0], bcts[:,1], label='bct')
181+
plt.scatter(ts[:,0], ts[:,1], label='sc')
182+
plt.scatter(ds[:,0], ds[:,1], label='diamond')
183+
#for j, a in enumerate(_as):
184+
# plt.plot(hcps[j,:,0], hcps[j,:,1], '--', label='hcp'+str(a))
185+
# plt.plot(bcts[j,:,0], bcts[j,:,1], '-d', label='bct'+str(a))
186+
# plt.plot(ts[j,:,0], ts[j,:,1], '-*', label='tetra'+str(a))
187+
# plt.plot(ds[j,:,0], ds[j,:,1], '-s', label='diamond'+str(a))
188+
plt.legend()
189+
plt.xlim([0,1])
190+
plt.ylim([0,1])
191+
plt.xlabel('q4')
192+
plt.ylabel('q6')
193+
plt.savefig('0.png')
302194

195+
303196

‎pyxtal/symmetry.py‎

Lines changed: 13 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -968,44 +968,18 @@ def from_group_and_index(group, index, dim=3, PBC=None):
968968
if c.isalpha():
969969
index = c
970970
break
971-
972-
if dim == 3:
973-
if wp.PBC is None:
974-
wp.PBC = [1, 1, 1]
975-
976-
ops_all = get_wyckoffs(wp.number)
977-
if use_letter:
978-
wp.index = index_from_letter(index, ops_all, dim=dim)
979-
wp.letter = index
980-
else:
981-
wp.letter = letter_from_index(wp.index, ops_all, dim=dim)
982-
wp.ops = ops_all[wp.index]
983-
984-
wp.multiplicity = len(wp.ops)
985-
wp.symmetry = get_wyckoff_symmetry(wp.number)[wp.index]
986-
wp.generators = get_generators(wp.number, dim)[wp.index]
987-
988-
elif dim == 2:
989-
if wp.PBC is None:
990-
wp.PBC = [1, 1, 0]
991-
992-
ops_all = get_layer(wp.number)
993-
if use_letter:
994-
wp.index = index_from_letter(index, ops_all, dim=dim)
995-
wp.letter = index
996-
else:
997-
wp.letter = letter_from_index(wp.index, ops_all, dim=dim)
998-
wp.ops = ops_all[wp.index]
999-
1000-
wp.multiplicity = len(wp.ops)
1001-
wp.symmetry = get_layer_symmetry(wp.number)[wp.index]
1002-
wp.generators = get_generators(wp.number, dim)[wp.index]
1003-
1004-
elif dim == 1:
1005-
if wp.PBC is None:
1006-
wp.PBC = [0, 0, 1]
1007-
1008-
ops_all = get_rod(wp.number)
971+
if dim > 0:
972+
if dim == 3:
973+
if wp.PBC is None:
974+
wp.PBC = [1, 1, 1]
975+
elif dim == 2:
976+
if wp.PBC is None:
977+
wp.PBC = [1, 1, 0]
978+
elif dim == 1:
979+
if wp.PBC is None:
980+
wp.PBC = [0, 0, 1]
981+
982+
ops_all = get_wyckoffs(wp.number, dim=dim)
1009983
if use_letter:
1010984
wp.index = index_from_letter(index, ops_all, dim=dim)
1011985
wp.letter = index
@@ -1014,7 +988,7 @@ def from_group_and_index(group, index, dim=3, PBC=None):
1014988
wp.ops = ops_all[wp.index]
1015989

1016990
wp.multiplicity = len(wp.ops)
1017-
wp.symmetry = get_rod_symmetry(wp.number)[wp.index]
991+
wp.symmetry = get_wyckoff_symmetry(wp.number, dim)[wp.index]
1018992
wp.generators = get_generators(wp.number, dim)[wp.index]
1019993

1020994
elif dim == 0:

0 commit comments

Comments
 (0)