1 | # -*- coding: utf-8 -*- |
---|
2 | ''' |
---|
3 | *GSASIIlattice: Unit cells* |
---|
4 | --------------------------- |
---|
5 | |
---|
6 | Perform lattice-related computations |
---|
7 | |
---|
8 | Note that *G* is the reciprocal lattice tensor, and *g* is its inverse, |
---|
9 | :math:`G = g^{-1}`, where |
---|
10 | |
---|
11 | .. math:: |
---|
12 | |
---|
13 | g = \\left( \\begin{matrix} |
---|
14 | a^2 & a b\\cos\\gamma & a c\\cos\\beta \\\\ |
---|
15 | a b\\cos\\gamma & b^2 & b c \\cos\\alpha \\\\ |
---|
16 | a c\\cos\\beta & b c \\cos\\alpha & c^2 |
---|
17 | \\end{matrix}\\right) |
---|
18 | |
---|
19 | The "*A* tensor" terms are defined as |
---|
20 | :math:`A = (\\begin{matrix} G_{11} & G_{22} & G_{33} & 2G_{12} & 2G_{13} & 2G_{23}\\end{matrix})` and *A* can be used in this fashion: |
---|
21 | :math:`d^* = \\sqrt {A_0 h^2 + A_1 k^2 + A_2 l^2 + A_3 hk + A_4 hl + A_5 kl}`, where |
---|
22 | *d* is the d-spacing, and :math:`d^*` is the reciprocal lattice spacing, |
---|
23 | :math:`Q = 2 \\pi d^* = 2 \\pi / d`. |
---|
24 | Note that GSAS-II variables ``p::Ai`` (i = 0, 1,... 5) and ``p`` is a phase number are |
---|
25 | used for the *Ai* values. See :func:`A2cell`, :func:`cell2A` for interconversion between A and |
---|
26 | unit cell parameters; :func:`cell2Gmat` :func:`Gmat2cell` for G and cell parameters. |
---|
27 | |
---|
28 | When the hydrostatic/elastic strain coefficients (*Dij*, :math:`D_{ij}`) are used, they are added to the |
---|
29 | *A* tensor terms (Ai, :math:`A_{i}`) so that A is redefined |
---|
30 | :math:`A = (\\begin{matrix} A_{0} + D_{11} & A_{1} + D_{22} & A_{2} + D_{33} & A_{3} + D_{12} & A_{4} + D_{13} & A_{5} + D_{23}\\end{matrix})`. See :func:`cellDijFill`. |
---|
31 | Note that GSAS-II variables ``p:h:Dij`` (i,j = 1, 2, 3) and ``p`` is a phase number |
---|
32 | and ``h`` a histogram number are used for the *Dij* values. |
---|
33 | ''' |
---|
34 | ########### SVN repository information ################### |
---|
35 | # $Date: 2023-01-17 13:31:50 +0000 (Tue, 17 Jan 2023) $ |
---|
36 | # $Author: toby $ |
---|
37 | # $Revision: 5478 $ |
---|
38 | # $URL: trunk/GSASIIlattice.py $ |
---|
39 | # $Id: GSASIIlattice.py 5478 2023-01-17 13:31:50Z toby $ |
---|
40 | ########### SVN repository information ################### |
---|
41 | from __future__ import division, print_function |
---|
42 | import math |
---|
43 | import time |
---|
44 | import copy |
---|
45 | import sys |
---|
46 | import random as ran |
---|
47 | import numpy as np |
---|
48 | import numpy.linalg as nl |
---|
49 | import scipy.special as spsp |
---|
50 | import GSASIIpath |
---|
51 | import GSASIImath as G2mth |
---|
52 | import GSASIIspc as G2spc |
---|
53 | import GSASIIElem as G2elem |
---|
54 | GSASIIpath.SetVersionNumber("$Revision: 5478 $") |
---|
55 | # trig functions in degrees |
---|
56 | sind = lambda x: np.sin(x*np.pi/180.) |
---|
57 | asind = lambda x: 180.*np.arcsin(x)/np.pi |
---|
58 | tand = lambda x: np.tan(x*np.pi/180.) |
---|
59 | atand = lambda x: 180.*np.arctan(x)/np.pi |
---|
60 | atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi |
---|
61 | cosd = lambda x: np.cos(x*np.pi/180.) |
---|
62 | acosd = lambda x: 180.*np.arccos(x)/np.pi |
---|
63 | rdsq2d = lambda x,p: round(1.0/np.sqrt(x),p) |
---|
64 | try: # fails on doc build |
---|
65 | rpd = np.pi/180. |
---|
66 | RSQ2PI = 1./np.sqrt(2.*np.pi) |
---|
67 | SQ2 = np.sqrt(2.) |
---|
68 | RSQPI = 1./np.sqrt(np.pi) |
---|
69 | R2pisq = 1./(2.*np.pi**2) |
---|
70 | except TypeError: |
---|
71 | pass |
---|
72 | nxs = np.newaxis |
---|
73 | |
---|
74 | def sec2HMS(sec): |
---|
75 | """Convert time in sec to H:M:S string |
---|
76 | |
---|
77 | :param sec: time in seconds |
---|
78 | :return: H:M:S string (to nearest 100th second) |
---|
79 | |
---|
80 | """ |
---|
81 | H = int(sec//3600) |
---|
82 | M = int(sec//60-H*60) |
---|
83 | S = sec-3600*H-60*M |
---|
84 | return '%d:%2d:%.2f'%(H,M,S) |
---|
85 | |
---|
86 | def rotdMat(angle,axis=0): |
---|
87 | """Prepare rotation matrix for angle in degrees about axis(=0,1,2) |
---|
88 | |
---|
89 | :param angle: angle in degrees |
---|
90 | :param axis: axis (0,1,2 = x,y,z) about which for the rotation |
---|
91 | :return: rotation matrix - 3x3 numpy array |
---|
92 | |
---|
93 | """ |
---|
94 | if axis == 2: |
---|
95 | return np.array([[cosd(angle),-sind(angle),0],[sind(angle),cosd(angle),0],[0,0,1]]) |
---|
96 | elif axis == 1: |
---|
97 | return np.array([[cosd(angle),0,-sind(angle)],[0,1,0],[sind(angle),0,cosd(angle)]]) |
---|
98 | else: |
---|
99 | return np.array([[1,0,0],[0,cosd(angle),-sind(angle)],[0,sind(angle),cosd(angle)]]) |
---|
100 | |
---|
101 | def rotdMat4(angle,axis=0): |
---|
102 | """Prepare rotation matrix for angle in degrees about axis(=0,1,2) with scaling for OpenGL |
---|
103 | |
---|
104 | :param angle: angle in degrees |
---|
105 | :param axis: axis (0,1,2 = x,y,z) about which for the rotation |
---|
106 | :return: rotation matrix - 4x4 numpy array (last row/column for openGL scaling) |
---|
107 | |
---|
108 | """ |
---|
109 | Mat = rotdMat(angle,axis) |
---|
110 | return np.concatenate((np.concatenate((Mat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0) |
---|
111 | |
---|
112 | def fillgmat(cell): |
---|
113 | """Compute lattice metric tensor from unit cell constants |
---|
114 | |
---|
115 | :param cell: tuple with a,b,c,alpha, beta, gamma (degrees) |
---|
116 | :return: 3x3 numpy array |
---|
117 | |
---|
118 | """ |
---|
119 | a,b,c,alp,bet,gam = cell |
---|
120 | g = np.array([ |
---|
121 | [a*a, a*b*cosd(gam), a*c*cosd(bet)], |
---|
122 | [a*b*cosd(gam), b*b, b*c*cosd(alp)], |
---|
123 | [a*c*cosd(bet) ,b*c*cosd(alp), c*c]]) |
---|
124 | return g |
---|
125 | |
---|
126 | def cell2Gmat(cell): |
---|
127 | """Compute real and reciprocal lattice metric tensor from unit cell constants |
---|
128 | |
---|
129 | :param cell: tuple with a,b,c,alpha, beta, gamma (degrees) |
---|
130 | :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays) |
---|
131 | |
---|
132 | """ |
---|
133 | g = fillgmat(cell) |
---|
134 | G = nl.inv(g) |
---|
135 | return G,g |
---|
136 | |
---|
137 | def A2Gmat(A,inverse=True): |
---|
138 | """Fill real & reciprocal metric tensor (G) from A. |
---|
139 | |
---|
140 | :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23] |
---|
141 | :param bool inverse: if True return both G and g; else just G |
---|
142 | :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays) |
---|
143 | |
---|
144 | """ |
---|
145 | G = np.array([ |
---|
146 | [A[0], A[3]/2., A[4]/2.], |
---|
147 | [A[3]/2.,A[1], A[5]/2.], |
---|
148 | [A[4]/2.,A[5]/2., A[2]]]) |
---|
149 | if inverse: |
---|
150 | g = nl.inv(G) |
---|
151 | return G,g |
---|
152 | else: |
---|
153 | return G |
---|
154 | |
---|
155 | def Gmat2A(G): |
---|
156 | """Extract A from reciprocal metric tensor (G) |
---|
157 | |
---|
158 | :param G: reciprocal maetric tensor (3x3 numpy array |
---|
159 | :return: A = [G11,G22,G33,2*G12,2*G13,2*G23] |
---|
160 | |
---|
161 | """ |
---|
162 | return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]] |
---|
163 | |
---|
164 | def cell2A(cell): |
---|
165 | """Obtain A = [G11,G22,G33,2*G12,2*G13,2*G23] from lattice parameters |
---|
166 | |
---|
167 | :param cell: [a,b,c,alpha,beta,gamma] (degrees) |
---|
168 | :return: G reciprocal metric tensor as 3x3 numpy array |
---|
169 | |
---|
170 | """ |
---|
171 | G,g = cell2Gmat(cell) |
---|
172 | return Gmat2A(G) |
---|
173 | |
---|
174 | def A2cell(A): |
---|
175 | """Compute unit cell constants from A |
---|
176 | |
---|
177 | :param A: [G11,G22,G33,2*G12,2*G13,2*G23] G - reciprocal metric tensor |
---|
178 | :return: a,b,c,alpha, beta, gamma (degrees) - lattice parameters |
---|
179 | |
---|
180 | """ |
---|
181 | G,g = A2Gmat(A) |
---|
182 | return Gmat2cell(g) |
---|
183 | |
---|
184 | def Gmat2cell(g): |
---|
185 | """Compute real/reciprocal lattice parameters from real/reciprocal metric tensor (g/G) |
---|
186 | The math works the same either way. |
---|
187 | |
---|
188 | :param g (or G): real (or reciprocal) metric tensor 3x3 array |
---|
189 | :return: a,b,c,alpha, beta, gamma (degrees) (or a*,b*,c*,alpha*,beta*,gamma* degrees) |
---|
190 | |
---|
191 | """ |
---|
192 | oldset = np.seterr('raise') |
---|
193 | a = np.sqrt(max(0,g[0][0])) |
---|
194 | b = np.sqrt(max(0,g[1][1])) |
---|
195 | c = np.sqrt(max(0,g[2][2])) |
---|
196 | alp = acosd(g[2][1]/(b*c)) |
---|
197 | bet = acosd(g[2][0]/(a*c)) |
---|
198 | gam = acosd(g[0][1]/(a*b)) |
---|
199 | np.seterr(**oldset) |
---|
200 | return a,b,c,alp,bet,gam |
---|
201 | |
---|
202 | def invcell2Gmat(invcell): |
---|
203 | """Compute real and reciprocal lattice metric tensor from reciprocal |
---|
204 | unit cell constants |
---|
205 | |
---|
206 | :param invcell: [a*,b*,c*,alpha*, beta*, gamma*] (degrees) |
---|
207 | :return: reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays) |
---|
208 | |
---|
209 | """ |
---|
210 | G = fillgmat(invcell) |
---|
211 | g = nl.inv(G) |
---|
212 | return G,g |
---|
213 | |
---|
214 | def cellDijFill(pfx,phfx,SGData,parmDict): |
---|
215 | '''Returns the filled-out reciprocal cell (A) terms |
---|
216 | from the parameter dictionaries corrected for Dij. |
---|
217 | |
---|
218 | :param str pfx: parameter prefix ("n::", where n is a phase number) |
---|
219 | :param dict SGdata: a symmetry object |
---|
220 | :param dict parmDict: a dictionary of parameters |
---|
221 | |
---|
222 | :returns: A,sigA where each is a list of six terms with the A terms |
---|
223 | ''' |
---|
224 | if pfx+'D11' not in parmDict: |
---|
225 | return None |
---|
226 | if SGData['SGLaue'] in ['-1',]: |
---|
227 | A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'], |
---|
228 | parmDict[pfx+'A2']+parmDict[phfx+'D33'], |
---|
229 | parmDict[pfx+'A3']+parmDict[phfx+'D12'],parmDict[pfx+'A4']+parmDict[phfx+'D13'], |
---|
230 | parmDict[pfx+'A5']+parmDict[phfx+'D23']] |
---|
231 | elif SGData['SGLaue'] in ['2/m',]: |
---|
232 | if SGData['SGUniq'] == 'a': |
---|
233 | A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'], |
---|
234 | parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,0,parmDict[pfx+'A5']+parmDict[phfx+'D23']] |
---|
235 | elif SGData['SGUniq'] == 'b': |
---|
236 | A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'], |
---|
237 | parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,parmDict[pfx+'A4']+parmDict[phfx+'D13'],0] |
---|
238 | else: |
---|
239 | A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'], |
---|
240 | parmDict[pfx+'A2']+parmDict[phfx+'D33'],parmDict[pfx+'A3']+parmDict[phfx+'D12'],0,0] |
---|
241 | elif SGData['SGLaue'] in ['mmm',]: |
---|
242 | A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'], |
---|
243 | parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,0,0] |
---|
244 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
---|
245 | A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'], |
---|
246 | parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,0,0] |
---|
247 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
---|
248 | A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'], |
---|
249 | parmDict[pfx+'A2']+parmDict[phfx+'D33'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],0,0] |
---|
250 | elif SGData['SGLaue'] in ['3R', '3mR']: |
---|
251 | A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'], |
---|
252 | parmDict[pfx+'A0']+parmDict[phfx+'D11'], |
---|
253 | parmDict[pfx+'A3']+parmDict[phfx+'D23'],parmDict[pfx+'A3']+parmDict[phfx+'D23'], |
---|
254 | parmDict[pfx+'A3']+parmDict[phfx+'D23']] |
---|
255 | elif SGData['SGLaue'] in ['m3m','m3']: |
---|
256 | A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'], |
---|
257 | parmDict[pfx+'A0']+parmDict[phfx+'D11'],0,0,0] |
---|
258 | return A |
---|
259 | |
---|
260 | def CellDijCorr(Cell,SGData,Data,hist): |
---|
261 | '''Returns the cell corrected for Dij values. |
---|
262 | |
---|
263 | :param list Cell: lattice parameters |
---|
264 | :param dict SGdata: a symmetry object |
---|
265 | :param dict Data: phase data structure; contains set of Dij values |
---|
266 | :param str hist: histogram name |
---|
267 | |
---|
268 | :returns: cell corrected for Dij values |
---|
269 | ''' |
---|
270 | A = cell2A(Cell) |
---|
271 | Dij = Data[hist]['HStrain'][0] |
---|
272 | newA = AplusDij(A,Dij,SGData) |
---|
273 | return A2cell(newA) |
---|
274 | |
---|
275 | def AplusDij(A,Dij,SGData): |
---|
276 | ''' returns the A corrected by Dij |
---|
277 | |
---|
278 | :param list A: reciprocal metric terms A0-A5 |
---|
279 | :param array Dij: unique Dij values as stored in Hstrain |
---|
280 | :param dict SGdata: a symmetry object |
---|
281 | |
---|
282 | :returns list newA: A corrected by Dij |
---|
283 | ''' |
---|
284 | if SGData['SGLaue'] in ['-1',]: |
---|
285 | newA = [A[0]+Dij[0],A[1]+Dij[1],A[2]+Dij[2],A[3]+Dij[3],A[4]+Dij[4],A[5]+Dij[5]] |
---|
286 | elif SGData['SGLaue'] in ['2/m',]: |
---|
287 | if SGData['SGUniq'] == 'a': |
---|
288 | newA = [A[0]+Dij[0],A[1]+Dij[1],A[2]+Dij[2],0,0,A[5]+Dij[3]] |
---|
289 | elif SGData['SGUniq'] == 'b': |
---|
290 | newA = [A[0]+Dij[0],A[1]+Dij[1],A[2]+Dij[2],0,A[4]+Dij[3],0] |
---|
291 | else: |
---|
292 | newA = [A[0]+Dij[0],A[1]+Dij[1],A[2]+Dij[2],A[3]+Dij[3],0,0] |
---|
293 | elif SGData['SGLaue'] in ['mmm',]: |
---|
294 | newA = [A[0]+Dij[0],A[1]+Dij[1],A[2]+Dij[2],0,0,0] |
---|
295 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
---|
296 | newA = [A[0]+Dij[0],A[0]+Dij[0],A[2]+Dij[1],0,0,0] |
---|
297 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
---|
298 | newA = [A[0]+Dij[0],A[0]+Dij[0],A[2]+Dij[1],A[0]+Dij[0],0,0] |
---|
299 | elif SGData['SGLaue'] in ['3R', '3mR']: |
---|
300 | newA = [A[0]+Dij[0],A[0]+Dij[0],A[0]+Dij[0],A[3]+Dij[1],A[3]+Dij[1],A[3]+Dij[1]] |
---|
301 | elif SGData['SGLaue'] in ['m3m','m3']: |
---|
302 | newA = [A[0]+Dij[0],A[0]+Dij[0],A[0]+Dij[0],0,0,0] |
---|
303 | |
---|
304 | return newA |
---|
305 | |
---|
306 | def prodMGMT(G,Mat): |
---|
307 | '''Transform metric tensor by matrix |
---|
308 | |
---|
309 | :param G: array metric tensor |
---|
310 | :param Mat: array transformation matrix |
---|
311 | :return: array new metric tensor |
---|
312 | |
---|
313 | ''' |
---|
314 | return np.inner(np.inner(Mat,G),Mat) #right |
---|
315 | # return np.inner(Mat,np.inner(Mat,G)) #right |
---|
316 | # return np.inner(np.inner(G,Mat).T,Mat) #right |
---|
317 | # return np.inner(Mat,np.inner(G,Mat).T) #right |
---|
318 | |
---|
319 | def TransformCell(cell,Trans): |
---|
320 | '''Transform lattice parameters by matrix |
---|
321 | |
---|
322 | :param cell: list a,b,c,alpha,beta,gamma,(volume) |
---|
323 | :param Trans: array transformation matrix |
---|
324 | :return: array transformed a,b,c,alpha,beta,gamma,volume |
---|
325 | |
---|
326 | ''' |
---|
327 | newCell = np.zeros(7) |
---|
328 | g = cell2Gmat(cell)[1] |
---|
329 | newg = prodMGMT(g,Trans) |
---|
330 | newCell[:6] = Gmat2cell(newg) |
---|
331 | newCell[6] = calc_V(cell2A(newCell[:6])) |
---|
332 | return newCell |
---|
333 | |
---|
334 | # code to generate lattice constraint relationships between two phases |
---|
335 | # (chemical & magnetic) related by a transformation matrix. |
---|
336 | |
---|
337 | def symInner(M1,M2): |
---|
338 | '''Compute inner product of two square matrices with symbolic processing |
---|
339 | Use dot product because sympy does not define an inner product primitive |
---|
340 | |
---|
341 | This requires that M1 & M2 be two sympy objects, as created in |
---|
342 | GenerateCellConstraints(). |
---|
343 | |
---|
344 | Note that this is only used to do the symbolic math needed to generate |
---|
345 | cell relationships. It is not used normally in GSAS-II. |
---|
346 | ''' |
---|
347 | import sympy as sym |
---|
348 | prodOuter = [] |
---|
349 | for i in range(3): |
---|
350 | prod = [] |
---|
351 | for j in range(3): |
---|
352 | prod.append(M1[i,:].dot(M2[j,:])) |
---|
353 | prodOuter.append(prod) |
---|
354 | return sym.Matrix(prodOuter) |
---|
355 | |
---|
356 | def GenerateCellConstraints(): |
---|
357 | '''Generate unit cell constraints for transforming one set of A tensor |
---|
358 | values to another using symbolic math (requires the sympy package) |
---|
359 | |
---|
360 | Note that this is only used to do the symbolic math needed to generate |
---|
361 | cell relationships. It is not used normally in GSAS-II. |
---|
362 | ''' |
---|
363 | import sympy as sym |
---|
364 | # define A tensor for starting cell |
---|
365 | A0, A1, A2, A3, A4, A5 = sym.symbols('A0, A1, A2, A3, A4, A5') |
---|
366 | G = sym.Matrix([ [A0, A3/2., A4/2.], |
---|
367 | [A3/2., A1, A5/2.], |
---|
368 | [A4/2., A5/2., A2 ]] ) |
---|
369 | # define transformation matrix |
---|
370 | T00, T10, T20, T01, T11, T21, T02, T12, T22 = sym.symbols( |
---|
371 | 'T00, T10, T20, T01, T11, T21, T02, T12, T22') |
---|
372 | Tr = sym.Matrix([ [T00, T10, T20], [T01, T11, T21], [T02, T12, T22],]) |
---|
373 | # apply transform |
---|
374 | newG = symInner(symInner(Tr,G),Tr).expand() |
---|
375 | # define A tensor for converted cell |
---|
376 | return [newG[0,0],newG[1,1],newG[2,2],2.*newG[0,1],2.*newG[0,2],2.*newG[1,2]] |
---|
377 | |
---|
378 | def subVals(expr,A,T): |
---|
379 | '''Evaluate the symbolic expressions by substituting for A0-A5 & Tij |
---|
380 | |
---|
381 | This can be used on the cell relationships created in |
---|
382 | :func:`GenerateCellConstraints` like this:: |
---|
383 | |
---|
384 | Trans = np.array([ [2/3, 4/3, 1/3], [-1, 0, 0], [-1/3, -2/3, 1/3] ]) |
---|
385 | T = np.linalg.inv(Trans).T |
---|
386 | print([subVals(i,Aold,T) for i in GenerateCellConstraints()]) |
---|
387 | |
---|
388 | :param list expr: a list of sympy expressions. |
---|
389 | :param list A: This is the A* tensor as defined above. |
---|
390 | :param np.array T: a 3x3 transformation matrix where, |
---|
391 | Trans = np.array([ [2/3, 4/3, 1/3], [-1, 0, 0], [-1/3, -2/3, 1/3] ]) |
---|
392 | (for a' = 2/3a + 4/3b + 1/3c; b' = -a; c' = -1/3, -2/3, 1/3) |
---|
393 | then T = np.linalg.inv(Trans).T |
---|
394 | |
---|
395 | Note that this is only used to do the symbolic math needed to generate |
---|
396 | cell relationships. It is not used normally in GSAS-II. |
---|
397 | ''' |
---|
398 | import sympy as sym |
---|
399 | A0, A1, A2, A3, A4, A5 = sym.symbols('A0, A1, A2, A3, A4, A5') |
---|
400 | # transformation matrix |
---|
401 | T00, T10, T20, T01, T11, T21, T02, T12, T22 = sym.symbols( |
---|
402 | 'T00, T10, T20, T01, T11, T21, T02, T12, T22') |
---|
403 | vals = dict(zip([T00, T10, T20, T01, T11, T21, T02, T12, T22],T.ravel())) |
---|
404 | vals.update(dict(zip([A0, A1, A2, A3, A4, A5],A))) |
---|
405 | return float(expr.subs(vals)) |
---|
406 | |
---|
407 | # some sample test code using the routines above follows:: |
---|
408 | # Trans = np.array([ [2/3, 4/3, 1/3], [-1, 0, 0], [-1/3, -2/3, 1/3] ]) |
---|
409 | # Mat = np.linalg.inv(Trans).T |
---|
410 | # Aold = [0.05259986634758891, 0.05259986634758891, 0.005290771904550856, 0.052599866347588925, 0, 0] |
---|
411 | # Anew = [0.018440738491448085, 0.03944989976069168, 0.034313054205100654, 0, -0.00513684555559103, 0] |
---|
412 | # cellConstr = G2lat.GenerateCellConstraints() |
---|
413 | # calcA = [G2lat.subVals(i,Aold,Mat) for i in cellConstr] |
---|
414 | # print('original xform A',Anew) |
---|
415 | # print('calculated xfrom A',calcA) |
---|
416 | # print('input') |
---|
417 | # print(' old cell',G2lat.A2cell(Aold)) |
---|
418 | # print(' new cell',G2lat.A2cell(Anew)) |
---|
419 | # print('derived results') |
---|
420 | # print(' from eq.',G2lat.A2cell(calcA)) |
---|
421 | # print(' diffs ',np.array(G2lat.A2cell(calcA)) - G2lat.A2cell(Anew)) |
---|
422 | |
---|
423 | def fmtCellConstraints(cellConstr): |
---|
424 | '''Format the cell relationships created in :func:`GenerateCellConstraints` |
---|
425 | in a format that can be used to generate constraints. |
---|
426 | |
---|
427 | Use:: |
---|
428 | |
---|
429 | cXforms = G2lat.fmtCellConstraints(G2lat.GenerateCellConstraints()) |
---|
430 | |
---|
431 | Note that this is only used to do the symbolic math needed to generate |
---|
432 | cell relationships. It is not used normally in GSAS-II. |
---|
433 | ''' |
---|
434 | import re |
---|
435 | import sympy as sym |
---|
436 | A3, A4, A5 = sym.symbols('A3, A4, A5') |
---|
437 | consDict = {} |
---|
438 | for num,cons in enumerate(cellConstr): |
---|
439 | cons = str(cons.factor(A3,A4,A5,deep=True).simplify()) |
---|
440 | cons = re.sub('T([0-2]?)([0-2]?)',r'T[\2,\1]',cons) # Tij to T[j,i] |
---|
441 | l = [] |
---|
442 | for i in str(cons).split('+'): |
---|
443 | if ')' in i: |
---|
444 | l[-1] += ' + ' + i.strip() |
---|
445 | else: |
---|
446 | l.append(i.strip()) |
---|
447 | print("\nA'{} = ".format(num),str(cons)) |
---|
448 | consDict[num] = l |
---|
449 | return consDict |
---|
450 | |
---|
451 | cellXformRelations = {0: ['1.0*A0*T[0,0]**2', |
---|
452 | '1.0*A1*T[0,1]**2', |
---|
453 | '1.0*A2*T[0,2]**2', |
---|
454 | '1.0*A3*T[0,0]*T[0,1]', |
---|
455 | '1.0*A4*T[0,0]*T[0,2]', |
---|
456 | '1.0*A5*T[0,1]*T[0,2]'], |
---|
457 | 1: ['1.0*A0*T[1,0]**2', |
---|
458 | '1.0*A1*T[1,1]**2', |
---|
459 | '1.0*A2*T[1,2]**2', |
---|
460 | '1.0*A3*T[1,0]*T[1,1]', |
---|
461 | '1.0*A4*T[1,0]*T[1,2]', |
---|
462 | '1.0*A5*T[1,1]*T[1,2]'], |
---|
463 | 2: ['1.0*A0*T[2,0]**2', |
---|
464 | '1.0*A1*T[2,1]**2', |
---|
465 | '1.0*A2*T[2,2]**2', |
---|
466 | '1.0*A3*T[2,0]*T[2,1]', |
---|
467 | '1.0*A4*T[2,0]*T[2,2]', |
---|
468 | '1.0*A5*T[2,1]*T[2,2]'], |
---|
469 | 3: ['2.0*A0*T[0,0]*T[1,0]', |
---|
470 | '2.0*A1*T[0,1]*T[1,1]', |
---|
471 | '2.0*A2*T[0,2]*T[1,2]', |
---|
472 | '1.0*A3*(T[0,0]*T[1,1] + T[1,0]*T[0,1])', |
---|
473 | '1.0*A4*(T[0,0]*T[1,2] + T[1,0]*T[0,2])', |
---|
474 | '1.0*A5*(T[0,1]*T[1,2] + T[1,1]*T[0,2])'], |
---|
475 | 4: ['2.0*A0*T[0,0]*T[2,0]', |
---|
476 | '2.0*A1*T[0,1]*T[2,1]', |
---|
477 | '2.0*A2*T[0,2]*T[2,2]', |
---|
478 | '1.0*A3*(T[0,0]*T[2,1] + T[2,0]*T[0,1])', |
---|
479 | '1.0*A4*(T[0,0]*T[2,2] + T[2,0]*T[0,2])', |
---|
480 | '1.0*A5*(T[0,1]*T[2,2] + T[2,1]*T[0,2])'], |
---|
481 | 5: ['2.0*A0*T[1,0]*T[2,0]', |
---|
482 | '2.0*A1*T[1,1]*T[2,1]', |
---|
483 | '2.0*A2*T[1,2]*T[2,2]', |
---|
484 | '1.0*A3*(T[1,0]*T[2,1] + T[2,0]*T[1,1])', |
---|
485 | '1.0*A4*(T[1,0]*T[2,2] + T[2,0]*T[1,2])', |
---|
486 | '1.0*A5*(T[1,1]*T[2,2] + T[2,1]*T[1,2])']} |
---|
487 | |
---|
488 | '''cellXformRelations provide the constraints on newA[i] values for a new |
---|
489 | cell generated from oldA[i] values. |
---|
490 | ''' |
---|
491 | # cellXformRelations values were generated using:: |
---|
492 | # from GSASIIlattice import fmtCellConstraints,GenerateCellConstraints |
---|
493 | # cellXformRelations = fmtCellConstraints(GenerateCellConstraints()) |
---|
494 | |
---|
495 | def GenCellConstraints(Trans,origPhase,newPhase,origA,oSGLaue,nSGLaue,debug=False): |
---|
496 | '''Generate the constraints between two unit cells constants for a phase transformed |
---|
497 | by matrix Trans. |
---|
498 | |
---|
499 | :param np.array Trans: a 3x3 direct cell transformation matrix where, |
---|
500 | Trans = np.array([ [2/3, 4/3, 1/3], [-1, 0, 0], [-1/3, -2/3, 1/3] ]) |
---|
501 | (for a' = 2/3a + 4/3b + 1/3c; b' = -a; c' = -1/3, -2/3, 1/3) |
---|
502 | :param int origPhase: phase id (pId) for original phase |
---|
503 | :param int newPhase: phase id for the transformed phase to be constrained from |
---|
504 | original phase |
---|
505 | :param list origA: reciprocal cell ("A*") tensor (used for debug only) |
---|
506 | :param dict oSGLaue: space group info for original phase |
---|
507 | :param dict nSGLaue: space group info for transformed phase |
---|
508 | :param bool debug: If true, the constraint input is used to compute and print A* |
---|
509 | and from that the direct cell for the transformed phase. |
---|
510 | ''' |
---|
511 | import GSASIIobj as G2obj |
---|
512 | T = Mat = np.linalg.inv(Trans).T |
---|
513 | Anew = [] |
---|
514 | constrList = [] |
---|
515 | uniqueAnew = cellUnique(nSGLaue) |
---|
516 | zeroAorig = cellZeros(oSGLaue) |
---|
517 | for i in range(6): |
---|
518 | constr = [[-1.0,G2obj.G2VarObj('{}::A{}'.format(newPhase,i))]] |
---|
519 | mult = [] |
---|
520 | for j,item in enumerate(cellXformRelations[i]): |
---|
521 | const, aTerm, tTerm = item.split('*',2) |
---|
522 | const = float(const) * eval(tTerm) |
---|
523 | mult.append(const) |
---|
524 | # skip over A terms that are required to be zero |
---|
525 | if zeroAorig[int(aTerm[1])]: continue # only add non-zero terms |
---|
526 | # ignore terms where either the Transform contribution is zero [= abs() < 1e-8] |
---|
527 | # If the multiplier term is zero I don't think this accidental |
---|
528 | # but since it will not change there is no reason to include that |
---|
529 | # term in any case |
---|
530 | if abs(const) < 1e-8: continue |
---|
531 | constr.append([const,G2obj.G2VarObj('{}::{}'.format(origPhase,aTerm))]) |
---|
532 | if i in uniqueAnew: |
---|
533 | constrList.append(constr + [0.0,None,'c']) |
---|
534 | if debug: Anew.append(np.dot(origA,mult)) |
---|
535 | if debug: |
---|
536 | print('xformed A* ',Anew) |
---|
537 | print('xformed cell',A2cell(Anew)) |
---|
538 | return constrList |
---|
539 | |
---|
540 | def cellUnique(SGData): |
---|
541 | '''Returns the indices for the unique A tensor terms |
---|
542 | based on the Laue class. |
---|
543 | Any terms that are determined from others or are zero are not included. |
---|
544 | |
---|
545 | :param dict SGdata: a symmetry object |
---|
546 | :returns: a list of 0 to 6 terms with indices of the unique A terms |
---|
547 | ''' |
---|
548 | if SGData['SGLaue'] in ['-1',]: |
---|
549 | return [0,1,2,3,4,5] |
---|
550 | elif SGData['SGLaue'] in ['2/m',]: |
---|
551 | if SGData['SGUniq'] == 'a': |
---|
552 | return [0,1,2,5] |
---|
553 | elif SGData['SGUniq'] == 'b': |
---|
554 | return [0,1,2,4] |
---|
555 | else: |
---|
556 | return [0,1,2,3] |
---|
557 | elif SGData['SGLaue'] in ['mmm',]: |
---|
558 | return [0,1,2] |
---|
559 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
---|
560 | return [0,2] |
---|
561 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
---|
562 | return [0,2] |
---|
563 | elif SGData['SGLaue'] in ['3R', '3mR']: |
---|
564 | return [0,3] |
---|
565 | elif SGData['SGLaue'] in ['m3m','m3']: |
---|
566 | return [0,] |
---|
567 | |
---|
568 | def cellZeros(SGData): |
---|
569 | '''Returns a list with the A terms required to be zero based on Laue symmetry |
---|
570 | |
---|
571 | :param dict SGdata: a symmetry object |
---|
572 | :returns: A list of six terms where the values are True if the |
---|
573 | A term must be zero, False otherwise. |
---|
574 | ''' |
---|
575 | if SGData['SGLaue'] in ['-1',]: |
---|
576 | return 6*[False] |
---|
577 | elif SGData['SGLaue'] in ['2/m',]: |
---|
578 | if SGData['SGUniq'] == 'a': |
---|
579 | return [False,False,False,True,True,False] |
---|
580 | elif SGData['SGUniq'] == 'b': |
---|
581 | return [False,False,False,True,False,True] |
---|
582 | else: |
---|
583 | return [False,False,False,False,True,True] |
---|
584 | elif SGData['SGLaue'] in ['mmm',]: |
---|
585 | return [False,False,False,True,True,True] |
---|
586 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
---|
587 | return [False,False,False,True,True,True] |
---|
588 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
---|
589 | return [False,False,False,False,True,True] |
---|
590 | elif SGData['SGLaue'] in ['3R', '3mR']: |
---|
591 | return 6*[False] |
---|
592 | elif SGData['SGLaue'] in ['m3m','m3']: |
---|
593 | return [False,False,False,True,True,True] |
---|
594 | |
---|
595 | def TransformXYZ(XYZ,Trans,Vec): |
---|
596 | return np.inner(XYZ,Trans)+Vec |
---|
597 | |
---|
598 | def TransformU6(U6,Trans): |
---|
599 | Uij = np.inner(Trans,np.inner(U6toUij(U6),Trans).T)/nl.det(Trans) |
---|
600 | return UijtoU6(Uij) |
---|
601 | |
---|
602 | def ExpandCell(Atoms,atCodes,cx,Trans): |
---|
603 | Unit = [int(max(abs(np.array(unit)))-1) for unit in Trans.T] |
---|
604 | nUnit = (Unit[0]+1)*(Unit[1]+1)*(Unit[2]+1) |
---|
605 | Ugrid = np.mgrid[0:Unit[0]+1,0:Unit[1]+1,0:Unit[2]+1] |
---|
606 | Ugrid = np.reshape(Ugrid,(3,nUnit)).T |
---|
607 | Codes = copy.deepcopy(atCodes) |
---|
608 | newAtoms = copy.deepcopy(Atoms) |
---|
609 | for unit in Ugrid[1:]: |
---|
610 | moreAtoms = copy.deepcopy(Atoms) |
---|
611 | for atom in moreAtoms: |
---|
612 | atom[cx:cx+3] += unit |
---|
613 | newAtoms += moreAtoms |
---|
614 | codes = copy.deepcopy(atCodes) |
---|
615 | moreCodes = [code+'+%d,%d,%d'%(unit[0],unit[1],unit[2]) for code in codes] |
---|
616 | Codes += moreCodes |
---|
617 | return newAtoms,Codes |
---|
618 | |
---|
619 | def TransformPhase(oldPhase,newPhase,Trans,Uvec,Vvec,ifMag,Force=True): |
---|
620 | '''Transform atoms from oldPhase to newPhase |
---|
621 | M' is inv(M) |
---|
622 | does X' = M(X-U)+V transformation for coordinates and U' = MUM/det(M) |
---|
623 | for anisotropic thermal parameters |
---|
624 | |
---|
625 | :param oldPhase: dict G2 phase info for old phase |
---|
626 | :param newPhase: dict G2 phase info for new phase; with new cell & space group |
---|
627 | atoms are from oldPhase & will be transformed |
---|
628 | :param Trans: lattice transformation matrix M |
---|
629 | :param Uvec: array parent coordinates transformation vector U |
---|
630 | :param Vvec: array child coordinate transformation vector V |
---|
631 | :param ifMag: bool True if convert to magnetic phase; |
---|
632 | if True all nonmagnetic atoms will be removed |
---|
633 | |
---|
634 | :return: newPhase dict modified G2 phase info |
---|
635 | :return: atCodes list atom transformation codes |
---|
636 | |
---|
637 | ''' |
---|
638 | cx,ct,cs,cia = oldPhase['General']['AtomPtrs'] |
---|
639 | cm = 0 |
---|
640 | if oldPhase['General']['Type'] == 'magnetic': |
---|
641 | cm = cx+4 |
---|
642 | oAmat,oBmat = cell2AB(oldPhase['General']['Cell'][1:7]) |
---|
643 | nAmat,nBmat = cell2AB(newPhase['General']['Cell'][1:7]) |
---|
644 | SGData = newPhase['General']['SGData'] |
---|
645 | invTrans = nl.inv(Trans) |
---|
646 | newAtoms,atCodes = FillUnitCell(oldPhase,Force) |
---|
647 | newAtoms,atCodes = ExpandCell(newAtoms,atCodes,cx,Trans) |
---|
648 | if ifMag: |
---|
649 | cia += 3 |
---|
650 | cs += 3 |
---|
651 | newPhase['General']['Type'] = 'magnetic' |
---|
652 | newPhase['General']['AtomPtrs'] = [cx,ct,cs,cia] |
---|
653 | magAtoms = [] |
---|
654 | magatCodes = [] |
---|
655 | Landeg = 2.0 |
---|
656 | for iat,atom in enumerate(newAtoms): |
---|
657 | if len(G2elem.GetMFtable([atom[ct],],[Landeg,])): |
---|
658 | magAtoms.append(atom[:cx+4]+[0.,0.,0.]+atom[cx+4:]) |
---|
659 | magatCodes.append(atCodes[iat]) |
---|
660 | newAtoms = magAtoms |
---|
661 | atCodes = magatCodes |
---|
662 | newPhase['Draw Atoms'] = [] |
---|
663 | for atom in newAtoms: |
---|
664 | xyz = TransformXYZ(atom[cx:cx+3]+Uvec,invTrans.T,Vvec) |
---|
665 | if Force: |
---|
666 | xyz = np.around(xyz,6)%1. |
---|
667 | atom[cx:cx+3] = xyz |
---|
668 | if atom[cia] == 'A': |
---|
669 | atom[cia+2:cia+8] = TransformU6(atom[cia+2:cia+8],Trans) |
---|
670 | atom[cs:cs+2] = G2spc.SytSym(atom[cx:cx+3],SGData)[:2] |
---|
671 | atom[cia+8] = ran.randint(0,sys.maxsize) |
---|
672 | if cm: |
---|
673 | mag = np.sqrt(np.sum(np.array(atom[cm:cm+3])**2)) |
---|
674 | if mag: |
---|
675 | mom = np.inner(np.array(atom[cm:cm+3]),oBmat) |
---|
676 | mom = np.inner(mom,invTrans) |
---|
677 | mom = np.inner(mom,nAmat) |
---|
678 | mom /= np.sqrt(np.sum(mom**2)) |
---|
679 | atom[cm:cm+3] = mom*mag |
---|
680 | newPhase['Atoms'] = newAtoms |
---|
681 | if SGData['SpGrp'] != 'P 1': |
---|
682 | newPhase['Atoms'],atCodes = GetUnique(newPhase,atCodes) |
---|
683 | newPhase['Drawing'] = [] |
---|
684 | newPhase['ranId'] = ran.randint(0,sys.maxsize) |
---|
685 | return newPhase,atCodes |
---|
686 | |
---|
687 | def FindNonstandard(controls,Phase): |
---|
688 | ''' |
---|
689 | Find nonstandard setting of magnetic cell that aligns with parent nuclear cell |
---|
690 | |
---|
691 | :param controls: list unit cell indexing controls |
---|
692 | :param Phase: dict new magnetic phase data (NB:not G2 phase construction); modified here |
---|
693 | :return: None |
---|
694 | |
---|
695 | ''' |
---|
696 | abc = np.eye(3) |
---|
697 | cba = np.rot90(np.eye(3)) |
---|
698 | cba[1,1] *= -1 #makes c-ba |
---|
699 | Mats = {'abc':abc,'cab':np.roll(abc,2,1),'bca':np.roll(abc,1,1), |
---|
700 | 'acb':np.roll(cba,1,1),'bac':np.roll(cba,2,1),'cba':cba} #ok |
---|
701 | BNS = {'A':{'abc':'A','cab':'C','bca':'B','acb':'A','bac':'B','cba':'C'}, |
---|
702 | 'B':{'abc':'B','cab':'A','bca':'C','acb':'C','bac':'A','cba':'B'}, |
---|
703 | 'C':{'abc':'C','cab':'B','bca':'A','acb':'B','bac':'C','cba':'A'}, |
---|
704 | 'a':{'abc':'a','cab':'c','bca':'b','acb':'a','bac':'b','cba':'c'}, #Ok |
---|
705 | 'b':{'abc':'b','cab':'a','bca':'c','acb':'c','bac':'a','cba':'b'}, |
---|
706 | 'c':{'abc':'c','cab':'b','bca':'a','acb':'b','bac':'c','cba':'a'}, |
---|
707 | 'S':{'abc':'S','cab':'S','bca':'S','acb':'S','bac':'S','cba':'S'}, |
---|
708 | 'I':{'abc':'I','cab':'I','bca':'I','acb':'I','bac':'I','cba':'I'}, |
---|
709 | } |
---|
710 | Trans = Phase['Trans'] |
---|
711 | Uvec = Phase['Uvec'] |
---|
712 | SGData = Phase['SGData'] |
---|
713 | MSG = SGData.get('MagSpGrp',SGData['SpGrp']).split(' ',1) |
---|
714 | MSG[0] += ' ' |
---|
715 | bns = '' |
---|
716 | if '_' in MSG[0]: |
---|
717 | bns = MSG[0][2] |
---|
718 | spn = SGData.get('SGSpin',[]) |
---|
719 | if 'ortho' in SGData['SGSys']: |
---|
720 | lattSym = G2spc.getlattSym(Trans) |
---|
721 | SpGrp = SGData['SpGrp'] |
---|
722 | NTrans = np.inner(Mats[lattSym].T,Trans.T) #ok |
---|
723 | if len(spn): spn[1:4] = np.inner(np.abs(nl.inv(Mats[lattSym])),spn[1:4]) #ok |
---|
724 | SGsym = G2spc.getlattSym(nl.inv(Mats[lattSym])) |
---|
725 | |
---|
726 | if lattSym != 'abc': |
---|
727 | NSG = G2spc.altSettingOrtho[SpGrp][SGsym].replace("'",'').split(' ') |
---|
728 | if ' '.join(NSG) in ['P 2 21 2',]: |
---|
729 | Uvec[1] += .25 |
---|
730 | elif ' '.join(NSG) in ['P 21 2 2',]: |
---|
731 | Uvec[0] += .25 |
---|
732 | elif ' '.join(NSG) in ['P 2 2 21',]: |
---|
733 | Uvec[2] += .25 |
---|
734 | Bns = '' |
---|
735 | if bns: |
---|
736 | Bns = BNS[bns][lattSym] |
---|
737 | NSG[0] += '_'+Bns+' ' |
---|
738 | elif len(spn): |
---|
739 | for ifld in [1,2,3]: |
---|
740 | if spn[ifld] < 0: |
---|
741 | NSG[ifld] += "'" |
---|
742 | Nresult = [''.join(NSG)+' ',Bns] |
---|
743 | return Nresult,Uvec,NTrans |
---|
744 | else: |
---|
745 | return None |
---|
746 | elif 'mono' in SGData['SGSys']: # and not 'P_A' in Phase['Name']: #skip the one that doesn't work |
---|
747 | newcell = TransformCell(controls[6:12],Trans) |
---|
748 | MatsA = np.array([[1.,0.,0.],[0.,1.,0.],[1.,0,1.]]) |
---|
749 | MatsB = np.array([[1.,0.,0.],[0.,1.,0.],[-1.,0,1.]]) |
---|
750 | if not 70. < newcell[4] < 110.: |
---|
751 | MSG[1] = MSG[1].replace('c','n') |
---|
752 | MSG[0] = MSG[0].replace('C_c','C_B').replace('P_A','P ') |
---|
753 | if '_' in MSG[0]: |
---|
754 | bns = MSG[0][2] |
---|
755 | if newcell[4] > 110.: |
---|
756 | if newcell[2] > newcell[0]: |
---|
757 | Mats = MatsA |
---|
758 | else: |
---|
759 | MSG[1] = MSG[1].replace('n','c') |
---|
760 | MSG[0] = MSG[0].replace('C ','I ') |
---|
761 | Mats = MatsA.T |
---|
762 | elif newcell[4] < 70.: |
---|
763 | if newcell[2] > newcell[0]: |
---|
764 | Mats = MatsB |
---|
765 | else: |
---|
766 | MSG[1] = MSG[1].replace('n','c') |
---|
767 | MSG[0] = MSG[0].replace('C ','I ') |
---|
768 | Mats = MatsB.T |
---|
769 | Nresult = [' '.join(MSG)+' ',bns] |
---|
770 | NTrans = np.inner(Mats,Trans.T) |
---|
771 | return Nresult,Uvec,NTrans |
---|
772 | return None |
---|
773 | |
---|
774 | def makeBilbaoPhase(result,uvec,trans,ifMag=False): |
---|
775 | phase = {} |
---|
776 | phase['Name'] = result[0].strip() |
---|
777 | phase['Uvec'] = uvec |
---|
778 | phase['Trans'] = trans |
---|
779 | phase['Keep'] = False |
---|
780 | phase['Use'] = False |
---|
781 | phase['aType'] = '' |
---|
782 | SpGp = result[0].replace("'",'') |
---|
783 | SpGrp = G2spc.StandardizeSpcName(SpGp) |
---|
784 | phase['SGData'] = G2spc.SpcGroup(SpGrp)[1] |
---|
785 | if ifMag: |
---|
786 | BNSlatt = phase['SGData']['SGLatt'] |
---|
787 | if not result[1]: |
---|
788 | phase['SGData']['SGSpin'] = G2spc.GetSGSpin(phase['SGData'],result[0]) |
---|
789 | phase['SGData']['GenSym'],phase['SGData']['GenFlg'],BNSsym = G2spc.GetGenSym(phase['SGData']) |
---|
790 | if result[1]: |
---|
791 | BNSlatt += '_'+result[1] |
---|
792 | if 'P_S' in BNSlatt: BNSlatt = 'P_c' #triclinic fix |
---|
793 | phase['SGData']['BNSlattsym'] = [BNSlatt,BNSsym[BNSlatt]] |
---|
794 | G2spc.ApplyBNSlatt(phase['SGData'],phase['SGData']['BNSlattsym']) |
---|
795 | phase['SGData']['SpnFlp'] = G2spc.GenMagOps(phase['SGData'])[1] |
---|
796 | phase['SGData']['MagSpGrp'] = G2spc.MagSGSym(phase['SGData']) |
---|
797 | return phase |
---|
798 | |
---|
799 | def FillUnitCell(Phase,Force=True): |
---|
800 | Atoms = copy.deepcopy(Phase['Atoms']) |
---|
801 | atomData = [] |
---|
802 | atCodes = [] |
---|
803 | SGData = Phase['General']['SGData'] |
---|
804 | SpnFlp = SGData.get('SpnFlp',[]) |
---|
805 | Amat,Bmat = cell2AB(Phase['General']['Cell'][1:7]) |
---|
806 | cx,ct,cs,cia = Phase['General']['AtomPtrs'] |
---|
807 | cm = 0 |
---|
808 | if Phase['General']['Type'] == 'magnetic': |
---|
809 | cm = cx+4 |
---|
810 | for iat,atom in enumerate(Atoms): |
---|
811 | XYZ = np.array(atom[cx:cx+3]) |
---|
812 | xyz = XYZ |
---|
813 | cellj = np.zeros(3,dtype=np.int32) |
---|
814 | if Force: |
---|
815 | xyz,cellj = G2spc.MoveToUnitCell(xyz) |
---|
816 | if atom[cia] == 'A': |
---|
817 | Uij = atom[cia+2:cia+8] |
---|
818 | result = G2spc.GenAtom(xyz,SGData,False,Uij,Force) |
---|
819 | for item in result: |
---|
820 | item = list(item) |
---|
821 | item[2] += cellj |
---|
822 | # if item[0][2] >= .95: item[0][2] -= 1. |
---|
823 | atom[cx:cx+3] = item[0] |
---|
824 | atom[cia+2:cia+8] = item[1] |
---|
825 | if cm: |
---|
826 | Opr = abs(item[2])%100 |
---|
827 | M = SGData['SGOps'][Opr-1][0] |
---|
828 | opNum = G2spc.GetOpNum(item[2],SGData) |
---|
829 | mom = np.inner(np.array(atom[cm:cm+3]),Bmat) |
---|
830 | atom[cm:cm+3] = np.inner(np.inner(mom,M),Amat)*nl.det(M)*SpnFlp[opNum-1] |
---|
831 | atCodes.append('%d:%s'%(iat,str(item[2]))) |
---|
832 | atomData.append(atom[:cia+9]) #not SS stuff |
---|
833 | else: |
---|
834 | result = G2spc.GenAtom(xyz,SGData,False,Move=Force) |
---|
835 | for item in result: |
---|
836 | item = list(item) |
---|
837 | item[2] += cellj |
---|
838 | # if item[0][2] >= .95: item[0][2] -= 1. |
---|
839 | atom[cx:cx+3] = item[0] |
---|
840 | if cm: |
---|
841 | Opr = abs(item[1])%100 |
---|
842 | M = SGData['SGOps'][Opr-1][0] |
---|
843 | opNum = G2spc.GetOpNum(item[1],SGData) |
---|
844 | mom = np.inner(np.array(atom[cm:cm+3]),Bmat) |
---|
845 | atom[cm:cm+3] = np.inner(np.inner(mom,M),Amat)*nl.det(M)*SpnFlp[opNum-1] |
---|
846 | atCodes.append('%d:%s'%(iat,str(item[1]))) |
---|
847 | atomData.append(atom[:cia+9]) #not SS stuff |
---|
848 | |
---|
849 | return atomData,atCodes |
---|
850 | |
---|
851 | def GetUnique(Phase,atCodes): |
---|
852 | |
---|
853 | def noDuplicate(xyzA,XYZ): |
---|
854 | if True in [np.allclose(xyzA%1.,xyzB%1.,atol=0.0002) for xyzB in XYZ]: |
---|
855 | return False |
---|
856 | return True |
---|
857 | |
---|
858 | cx,ct = Phase['General']['AtomPtrs'][:2] |
---|
859 | SGData = Phase['General']['SGData'] |
---|
860 | Atoms = Phase['Atoms'] |
---|
861 | Ind = len(Atoms) |
---|
862 | newAtoms = [] |
---|
863 | newAtCodes = [] |
---|
864 | Indx = {} |
---|
865 | XYZ = {} |
---|
866 | for ind in range(Ind): |
---|
867 | XYZ[ind] = np.array(Atoms[ind][cx:cx+3])%1. |
---|
868 | Indx[ind] = True |
---|
869 | for ind in range(Ind): |
---|
870 | if Indx[ind]: |
---|
871 | xyz = XYZ[ind] |
---|
872 | for jnd in range(Ind): |
---|
873 | if Atoms[ind][ct-1] == Atoms[jnd][ct-1]: |
---|
874 | if ind != jnd and Indx[jnd]: |
---|
875 | Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True) |
---|
876 | xyzs = np.array([equiv[0] for equiv in Equiv]) |
---|
877 | Indx[jnd] = noDuplicate(xyz,xyzs) |
---|
878 | Ind = [] |
---|
879 | for ind in Indx: |
---|
880 | if Indx[ind]: |
---|
881 | newAtoms.append(Atoms[ind]) |
---|
882 | newAtCodes.append(atCodes[ind]) |
---|
883 | return newAtoms,newAtCodes |
---|
884 | |
---|
885 | def calc_rVsq(A): |
---|
886 | """Compute the square of the reciprocal lattice volume (1/V**2) from A' |
---|
887 | |
---|
888 | """ |
---|
889 | G,g = A2Gmat(A) |
---|
890 | rVsq = nl.det(G) |
---|
891 | if rVsq < 0: |
---|
892 | return 1 |
---|
893 | return rVsq |
---|
894 | |
---|
895 | def calc_rV(A): |
---|
896 | """Compute the reciprocal lattice volume (V*) from A |
---|
897 | """ |
---|
898 | return np.sqrt(calc_rVsq(A)) |
---|
899 | |
---|
900 | def calc_V(A): |
---|
901 | """Compute the real lattice volume (V) from A |
---|
902 | """ |
---|
903 | return 1./calc_rV(A) |
---|
904 | |
---|
905 | def A2invcell(A): |
---|
906 | """Compute reciprocal unit cell constants from A |
---|
907 | returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees) |
---|
908 | """ |
---|
909 | G,g = A2Gmat(A) |
---|
910 | return Gmat2cell(G) |
---|
911 | |
---|
912 | def Gmat2AB(G): |
---|
913 | """Computes orthogonalization matrix from reciprocal metric tensor G |
---|
914 | |
---|
915 | :returns: tuple of two 3x3 numpy arrays (A,B) |
---|
916 | |
---|
917 | * A for crystal to Cartesian transformations (A*x = np.inner(A,x) = X) |
---|
918 | * B (= inverse of A) for Cartesian to crystal transformation (B*X = np.inner(B,X) = x) |
---|
919 | |
---|
920 | """ |
---|
921 | # cellstar = Gmat2cell(G) |
---|
922 | g = nl.inv(G) |
---|
923 | cell = Gmat2cell(g) |
---|
924 | # A = np.zeros(shape=(3,3)) |
---|
925 | return cell2AB(cell) |
---|
926 | # # from Giacovazzo (Fundamentals 2nd Ed.) p.75 |
---|
927 | # A[0][0] = cell[0] # a |
---|
928 | # A[0][1] = cell[1]*cosd(cell[5]) # b cos(gamma) |
---|
929 | # A[0][2] = cell[2]*cosd(cell[4]) # c cos(beta) |
---|
930 | # A[1][1] = cell[1]*sind(cell[5]) # b sin(gamma) |
---|
931 | # A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta) |
---|
932 | # A[2][2] = 1./cellstar[2] # 1/c* |
---|
933 | # B = nl.inv(A) |
---|
934 | # return A,B |
---|
935 | |
---|
936 | def cell2AB(cell,alt=False): |
---|
937 | """Computes orthogonalization matrix from unit cell constants |
---|
938 | |
---|
939 | :param tuple cell: a,b,c, alpha, beta, gamma (degrees) |
---|
940 | :returns: tuple of two 3x3 numpy arrays (A,B) |
---|
941 | A for crystal to Cartesian transformations A*x = np.inner(A,x) = X |
---|
942 | B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B,X) = x |
---|
943 | """ |
---|
944 | G,g = cell2Gmat(cell) |
---|
945 | cellstar = Gmat2cell(G) |
---|
946 | A = np.zeros(shape=(3,3)) |
---|
947 | if alt: #as used in RMCProfile!! |
---|
948 | A[0][0] = 1./cellstar[0] |
---|
949 | A[0][1] = cell[0]*cosd(cell[5])*sind(cell[3]) |
---|
950 | A[0][2] = cell[0]*cosd(cell[4]) |
---|
951 | A[1][1] = cell[1]*sind(cell[3]) |
---|
952 | A[1][2] = cell[1]*cosd(cell[3]) |
---|
953 | A[2][2] = cell[2] |
---|
954 | B = nl.inv(A) |
---|
955 | return A,B |
---|
956 | # from Giacovazzo (Fundamentals 2nd Ed.) p.75 |
---|
957 | A[0][0] = cell[0] # a |
---|
958 | A[0][1] = cell[1]*cosd(cell[5]) # b cos(gamma) |
---|
959 | A[0][2] = cell[2]*cosd(cell[4]) # c cos(beta) |
---|
960 | A[1][1] = cell[1]*sind(cell[5]) # b sin(gamma) |
---|
961 | A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta) |
---|
962 | A[2][2] = 1./cellstar[2] # 1/c* |
---|
963 | B = nl.inv(A) |
---|
964 | return A,B |
---|
965 | |
---|
966 | def HKL2SpAng(H,cell,SGData): |
---|
967 | """Computes spherical coords for hkls; view along 001 |
---|
968 | |
---|
969 | :param array H: arrays of hkl |
---|
970 | :param tuple cell: a,b,c, alpha, beta, gamma (degrees) |
---|
971 | :param dict SGData: space group dictionary |
---|
972 | :returns: arrays of r,phi,psi (radius,inclination,azimuth) about 001 |
---|
973 | """ |
---|
974 | A,B = cell2AB(cell) |
---|
975 | xH = np.inner(B.T,H) |
---|
976 | r = np.sqrt(np.sum(xH**2,axis=0)) |
---|
977 | phi = acosd(xH[2]/r) |
---|
978 | psi = atan2d(xH[1],xH[0]) |
---|
979 | phi = np.where(phi>90.,180.-phi,phi) |
---|
980 | # GSASIIpath.IPyBreak() |
---|
981 | return r,phi,psi |
---|
982 | |
---|
983 | def U6toUij(U6): |
---|
984 | """Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23] |
---|
985 | NB: there is a non numpy version in GSASIIspc: U2Uij |
---|
986 | |
---|
987 | :param list U6: 6 terms of u11,u22,... |
---|
988 | :returns: |
---|
989 | Uij - numpy [3][3] array of uij |
---|
990 | """ |
---|
991 | U = np.array([ |
---|
992 | [U6[0], U6[3], U6[4]], |
---|
993 | [U6[3], U6[1], U6[5]], |
---|
994 | [U6[4], U6[5], U6[2]]]) |
---|
995 | return U |
---|
996 | |
---|
997 | def UijtoU6(U): |
---|
998 | """Fill vector [U11,U22,U33,U12,U13,U23] from Uij |
---|
999 | NB: there is a non numpy version in GSASIIspc: Uij2U |
---|
1000 | """ |
---|
1001 | U6 = np.array([U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]) |
---|
1002 | return U6 |
---|
1003 | |
---|
1004 | def betaij2Uij(betaij,G): |
---|
1005 | """ |
---|
1006 | Convert beta-ij to Uij tensors |
---|
1007 | |
---|
1008 | :param beta-ij - numpy array [beta-ij] |
---|
1009 | :param G: reciprocal metric tensor |
---|
1010 | :returns: Uij: numpy array [Uij] |
---|
1011 | """ |
---|
1012 | ast = np.sqrt(np.diag(G)) #a*, b*, c* |
---|
1013 | Mast = np.multiply.outer(ast,ast) |
---|
1014 | return R2pisq*UijtoU6(U6toUij(betaij)/Mast) |
---|
1015 | |
---|
1016 | def Uij2betaij(Uij,G): |
---|
1017 | """ |
---|
1018 | Convert Uij to beta-ij tensors -- stub for eventual completion |
---|
1019 | |
---|
1020 | :param Uij: numpy array [Uij] |
---|
1021 | :param G: reciprocal metric tensor |
---|
1022 | :returns: beta-ij - numpy array [beta-ij] |
---|
1023 | """ |
---|
1024 | pass |
---|
1025 | |
---|
1026 | def cell2GS(cell): |
---|
1027 | ''' returns Uij to betaij conversion matrix''' |
---|
1028 | G,g = cell2Gmat(cell) |
---|
1029 | GS = G |
---|
1030 | GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1]) |
---|
1031 | GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2]) |
---|
1032 | GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2]) |
---|
1033 | return GS |
---|
1034 | |
---|
1035 | def Uij2Ueqv(Uij,GS,Amat): |
---|
1036 | ''' returns 1/3 trace of diagonalized U matrix |
---|
1037 | :param Uij: numpy array [Uij] |
---|
1038 | :param GS: Uij too betaij conversion matrix |
---|
1039 | :param Amat: crystal to Cartesian transformation matrix |
---|
1040 | :returns: 1/3 trace of diagonalized U matrix |
---|
1041 | :returns: True if nonpositive-definate; False otherwise |
---|
1042 | ''' |
---|
1043 | U = np.multiply(U6toUij(Uij),GS) |
---|
1044 | U = np.inner(Amat,np.inner(U,Amat).T) |
---|
1045 | E,R = nl.eigh(U) |
---|
1046 | return np.sum(E)/3.,E[0] < 0. |
---|
1047 | |
---|
1048 | def CosAngle(U,V,G): |
---|
1049 | """ calculate cos of angle between U & V in generalized coordinates |
---|
1050 | defined by metric tensor G |
---|
1051 | |
---|
1052 | :param U: 3-vectors assume numpy arrays, can be multiple reflections as (N,3) array |
---|
1053 | :param V: 3-vectors assume numpy arrays, only as (3) vector |
---|
1054 | :param G: metric tensor for U & V defined space assume numpy array |
---|
1055 | :returns: |
---|
1056 | cos(phi) |
---|
1057 | """ |
---|
1058 | u = (U.T/np.sqrt(np.sum(np.inner(U,G)*U,axis=1))).T |
---|
1059 | v = V/np.sqrt(np.inner(V,np.inner(G,V))) |
---|
1060 | cosP = np.inner(u,np.inner(G,v)) |
---|
1061 | return cosP |
---|
1062 | |
---|
1063 | def CosSinAngle(U,V,G): |
---|
1064 | """ calculate sin & cos of angle between U & V in generalized coordinates |
---|
1065 | defined by metric tensor G |
---|
1066 | |
---|
1067 | :param U: 3-vectors assume numpy arrays |
---|
1068 | :param V: 3-vectors assume numpy arrays |
---|
1069 | :param G: metric tensor for U & V defined space assume numpy array |
---|
1070 | :returns: |
---|
1071 | cos(phi) & sin(phi) |
---|
1072 | """ |
---|
1073 | u = U/np.sqrt(np.inner(U,np.inner(G,U))) |
---|
1074 | v = V/np.sqrt(np.inner(V,np.inner(G,V))) |
---|
1075 | cosP = np.inner(u,np.inner(G,v)) |
---|
1076 | sinP = np.sqrt(max(0.0,1.0-cosP**2)) |
---|
1077 | return cosP,sinP |
---|
1078 | |
---|
1079 | def criticalEllipse(prob): |
---|
1080 | """ |
---|
1081 | Calculate critical values for probability ellipsoids from probability |
---|
1082 | """ |
---|
1083 | if not ( 0.01 <= prob < 1.0): |
---|
1084 | return 1.54 |
---|
1085 | coeff = np.array([6.44988E-09,4.16479E-07,1.11172E-05,1.58767E-04,0.00130554, |
---|
1086 | 0.00604091,0.0114921,-0.040301,-0.6337203,1.311582]) |
---|
1087 | llpr = math.log(-math.log(prob)) |
---|
1088 | return np.polyval(coeff,llpr) |
---|
1089 | |
---|
1090 | def CellBlock(nCells): |
---|
1091 | """ |
---|
1092 | Generate block of unit cells n*n*n on a side; [0,0,0] centered, n = 2*nCells+1 |
---|
1093 | currently only works for nCells = 0 or 1 (not >1) |
---|
1094 | """ |
---|
1095 | if nCells: |
---|
1096 | N = 2*nCells+1 |
---|
1097 | N2 = N*N |
---|
1098 | N3 = N*N*N |
---|
1099 | cellArray = [] |
---|
1100 | A = np.array(range(N3)) |
---|
1101 | cellGen = np.array([A//N2-1,A//N%N-1,A%N-1]).T |
---|
1102 | for cell in cellGen: |
---|
1103 | cellArray.append(cell) |
---|
1104 | return cellArray |
---|
1105 | else: |
---|
1106 | return [0,0,0] |
---|
1107 | |
---|
1108 | def CellAbsorption(ElList,Volume): |
---|
1109 | '''Compute unit cell absorption |
---|
1110 | |
---|
1111 | :param dict ElList: dictionary of element contents including mu and |
---|
1112 | number of atoms be cell |
---|
1113 | :param float Volume: unit cell volume |
---|
1114 | :returns: mu-total/Volume |
---|
1115 | ''' |
---|
1116 | muT = 0 |
---|
1117 | for El in ElList: |
---|
1118 | muT += ElList[El]['mu']*ElList[El]['FormulaNo'] |
---|
1119 | return muT/Volume |
---|
1120 | |
---|
1121 | #Permutations and Combinations |
---|
1122 | # Four routines: combinations,uniqueCombinations, selections & permutations |
---|
1123 | #These taken from Python Cookbook, 2nd Edition. 19.15 p724-726 |
---|
1124 | # |
---|
1125 | def _combinators(_handle, items, n): |
---|
1126 | """ factored-out common structure of all following combinators """ |
---|
1127 | if n==0: |
---|
1128 | yield [ ] |
---|
1129 | return |
---|
1130 | for i, item in enumerate(items): |
---|
1131 | this_one = [ item ] |
---|
1132 | for cc in _combinators(_handle, _handle(items, i), n-1): |
---|
1133 | yield this_one + cc |
---|
1134 | def combinations(items, n): |
---|
1135 | """ take n distinct items, order matters """ |
---|
1136 | def skipIthItem(items, i): |
---|
1137 | return items[:i] + items[i+1:] |
---|
1138 | return _combinators(skipIthItem, items, n) |
---|
1139 | def uniqueCombinations(items, n): |
---|
1140 | """ take n distinct items, order is irrelevant """ |
---|
1141 | def afterIthItem(items, i): |
---|
1142 | return items[i+1:] |
---|
1143 | return _combinators(afterIthItem, items, n) |
---|
1144 | def selections(items, n): |
---|
1145 | """ take n (not necessarily distinct) items, order matters """ |
---|
1146 | def keepAllItems(items, i): |
---|
1147 | return items |
---|
1148 | return _combinators(keepAllItems, items, n) |
---|
1149 | def permutations(items): |
---|
1150 | """ take all items, order matters """ |
---|
1151 | return combinations(items, len(items)) |
---|
1152 | |
---|
1153 | #reflection generation routines |
---|
1154 | #for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor; |
---|
1155 | # cell - a,b,c,alp,bet,gam in A & deg |
---|
1156 | |
---|
1157 | def Pos2dsp(Inst,pos): |
---|
1158 | ''' convert powder pattern position (2-theta or TOF, musec) to d-spacing |
---|
1159 | is currently only approximate for EDX data; accurate for others. |
---|
1160 | ''' |
---|
1161 | if 'T' in Inst['Type'][0]: |
---|
1162 | return TOF2dsp(Inst,pos) |
---|
1163 | elif 'E' in Inst['Type'][0]: |
---|
1164 | return 12.398/(2.0*pos*sind(Inst['2-theta'][1]/2.0)) |
---|
1165 | else: #'PKS', 'C' or 'B' |
---|
1166 | wave = G2mth.getWave(Inst) |
---|
1167 | return wave/(2.0*sind((pos-Inst.get('Zero',[0,0])[1])/2.0)) |
---|
1168 | |
---|
1169 | def TOF2dsp(Inst,Pos): |
---|
1170 | ''' convert powder pattern TOF, musec to d-spacing by successive approximation |
---|
1171 | Pos can be numpy array |
---|
1172 | ''' |
---|
1173 | def func(d,pos,Inst): |
---|
1174 | return (pos-Inst['difA'][1]*d**2-Inst['Zero'][1]-Inst['difB'][1]/d)/Inst['difC'][1] |
---|
1175 | dsp0 = Pos/Inst['difC'][1] |
---|
1176 | N = 0 |
---|
1177 | while True: #successive approximations |
---|
1178 | dsp = func(dsp0,Pos,Inst) |
---|
1179 | if np.allclose(dsp,dsp0,atol=0.000001): |
---|
1180 | return dsp |
---|
1181 | dsp0 = dsp |
---|
1182 | N += 1 |
---|
1183 | if N > 10: |
---|
1184 | return dsp |
---|
1185 | |
---|
1186 | def Dsp2pos(Inst,dsp): |
---|
1187 | ''' convert d-spacing to powder pattern position (2-theta or TOF, musec) |
---|
1188 | ''' |
---|
1189 | if 'T' in Inst['Type'][0]: |
---|
1190 | pos = Inst['difC'][1]*dsp+Inst['Zero'][1]+Inst['difA'][1]*dsp**2+Inst.get('difB',[0,0,False])[1]/dsp |
---|
1191 | elif 'E' in Inst['Type'][0]: |
---|
1192 | return 12.398/(2.0*dsp*sind(Inst['2-theta'][1]/2.0))+Inst['ZE'][1]+Inst['YE'][1]*dsp+Inst['XE'][1]*dsp**2 |
---|
1193 | else: #'C' or 'B' |
---|
1194 | wave = G2mth.getWave(Inst) |
---|
1195 | val = min(0.995,wave/(2.*dsp)) #set max at 168deg |
---|
1196 | pos = 2.0*asind(val)+Inst.get('Zero',[0,0])[1] |
---|
1197 | return pos |
---|
1198 | |
---|
1199 | def getPeakPos(dataType,parmdict,dsp): |
---|
1200 | ''' convert d-spacing to powder pattern position (2-theta, E or TOF, musec) |
---|
1201 | ''' |
---|
1202 | if 'T' in dataType: |
---|
1203 | pos = parmdict['difC']*dsp+parmdict['difA']*dsp**2+parmdict['difB']/dsp+parmdict['Zero'] |
---|
1204 | elif 'E'in dataType: |
---|
1205 | pos = 12.398/(2.0*dsp*sind(parmdict['2-theta']/2.0)+parmdict['ZE']+parmdict['YE']*dsp+parmdict['XE']*dsp**2) |
---|
1206 | else: #'C' or 'B' |
---|
1207 | pos = 2.0*asind(parmdict['Lam']/(2.*dsp))+parmdict['Zero'] |
---|
1208 | return pos |
---|
1209 | |
---|
1210 | def calc_rDsq(H,A): |
---|
1211 | 'needs doc string' |
---|
1212 | rdsq = H[0]*H[0]*A[0]+H[1]*H[1]*A[1]+H[2]*H[2]*A[2]+H[0]*H[1]*A[3]+H[0]*H[2]*A[4]+H[1]*H[2]*A[5] |
---|
1213 | return rdsq |
---|
1214 | |
---|
1215 | def calc_rDsq2(H,G): |
---|
1216 | 'needs doc string' |
---|
1217 | return np.inner(H,np.inner(G,H)) |
---|
1218 | |
---|
1219 | def calc_rDsqSS(H,A,vec): |
---|
1220 | 'needs doc string' |
---|
1221 | rdsq = calc_rDsq(H[:3]+(H[3]*vec).T,A) |
---|
1222 | return rdsq |
---|
1223 | |
---|
1224 | def calc_rDsqZ(H,A,Z,tth,lam): |
---|
1225 | 'needs doc string' |
---|
1226 | rdsq = calc_rDsq(H,A)+Z*sind(tth)*2.0*rpd/lam**2 |
---|
1227 | return rdsq |
---|
1228 | |
---|
1229 | def calc_rDsqZSS(H,A,vec,Z,tth,lam): |
---|
1230 | 'needs doc string' |
---|
1231 | rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]*vec).T,A)+Z*sind(tth)*2.0*rpd/lam**2 |
---|
1232 | return rdsq |
---|
1233 | |
---|
1234 | def calc_rDsqT(H,A,Z,tof,difC): |
---|
1235 | 'needs doc string' |
---|
1236 | rdsq = calc_rDsq(H,A)+Z/difC |
---|
1237 | return rdsq |
---|
1238 | |
---|
1239 | def calc_rDsqTSS(H,A,vec,Z,tof,difC): |
---|
1240 | 'needs doc string' |
---|
1241 | rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]*vec).T,A)+Z/difC |
---|
1242 | return rdsq |
---|
1243 | |
---|
1244 | def PlaneIntercepts(Amat,H,phase,stack): |
---|
1245 | ''' find unit cell intercepts for a stack of hkl planes |
---|
1246 | ''' |
---|
1247 | Steps = range(-1,2,2) |
---|
1248 | if stack: |
---|
1249 | Steps = range(-10,10,1) |
---|
1250 | Stack = [] |
---|
1251 | Ux = np.array([[0,0],[1,0],[1,1],[0,1]]) |
---|
1252 | for step in Steps: |
---|
1253 | HX = [] |
---|
1254 | for i in [0,1,2]: |
---|
1255 | if H[i]: |
---|
1256 | h,k,l = [(i+1)%3,(i+2)%3,(i+3)%3] |
---|
1257 | for j in [0,1,2,3]: |
---|
1258 | hx = [0,0,0] |
---|
1259 | intcpt = ((phase)/360.+step-H[h]*Ux[j,0]-H[k]*Ux[j,1])/H[l] |
---|
1260 | if 0. <= intcpt <= 1.: |
---|
1261 | hx[h] = Ux[j,0] |
---|
1262 | hx[k] = Ux[j,1] |
---|
1263 | hx[l] = intcpt |
---|
1264 | HX.append(hx) |
---|
1265 | if len(HX)> 2: |
---|
1266 | HX = np.array(HX) |
---|
1267 | DX = np.inner(HX-HX[0],Amat) |
---|
1268 | D = np.sqrt(np.sum(DX**2,axis=1)) |
---|
1269 | Dsort = np.argsort(D) |
---|
1270 | HX = HX[Dsort] |
---|
1271 | DX = DX[Dsort] |
---|
1272 | D = D[Dsort] |
---|
1273 | DX[1:,:] = DX[1:,:]/D[1:,nxs] |
---|
1274 | A = 2.*np.ones(HX.shape[0]) |
---|
1275 | A[1:] = [np.dot(DX[1],dx) for dx in DX[1:]] |
---|
1276 | HX = HX[np.argsort(A)] |
---|
1277 | Stack.append(HX) |
---|
1278 | return Stack |
---|
1279 | |
---|
1280 | def MaxIndex(dmin,A): |
---|
1281 | 'needs doc string' |
---|
1282 | Hmax = [0,0,0] |
---|
1283 | try: |
---|
1284 | cell = A2cell(A) |
---|
1285 | except: |
---|
1286 | cell = [1.,1.,1.,90.,90.,90.] |
---|
1287 | for i in range(3): |
---|
1288 | Hmax[i] = int(np.round(cell[i]/dmin)) |
---|
1289 | return Hmax |
---|
1290 | |
---|
1291 | def transposeHKLF(transMat,Super,refList): |
---|
1292 | ''' Apply transformation matrix to hkl(m) |
---|
1293 | param: transmat: 3x3 or 4x4 array |
---|
1294 | param: Super: 0 or 1 for extra index |
---|
1295 | param: refList list of h,k,l,.... |
---|
1296 | return: newRefs transformed list of h',k',l',,, |
---|
1297 | return: badRefs list of noninteger h',k',l'... |
---|
1298 | ''' |
---|
1299 | newRefs = np.copy(refList) |
---|
1300 | badRefs = [] |
---|
1301 | for H in newRefs: |
---|
1302 | newH = np.inner(transMat,H[:3+Super]) |
---|
1303 | H[:3+Super] = np.rint(newH) |
---|
1304 | if not np.allclose(newH,H[:3+Super],atol=0.01): |
---|
1305 | badRefs.append(newH) |
---|
1306 | return newRefs,badRefs |
---|
1307 | |
---|
1308 | def sortHKLd(HKLd,ifreverse,ifdup,ifSS=False): |
---|
1309 | '''sort reflection list on d-spacing; can sort in either order |
---|
1310 | |
---|
1311 | :param HKLd: a list of [h,k,l,d,...]; |
---|
1312 | :param ifreverse: True for largest d first |
---|
1313 | :param ifdup: True if duplicate d-spacings allowed |
---|
1314 | :return: sorted reflection list |
---|
1315 | ''' |
---|
1316 | T = [] |
---|
1317 | N = 3 |
---|
1318 | if ifSS: |
---|
1319 | N = 4 |
---|
1320 | for i,H in enumerate(HKLd): |
---|
1321 | if ifdup: |
---|
1322 | T.append((H[N],i)) |
---|
1323 | else: |
---|
1324 | T.append(H[N]) |
---|
1325 | D = dict(zip(T,HKLd)) |
---|
1326 | T.sort() |
---|
1327 | if ifreverse: |
---|
1328 | T.reverse() |
---|
1329 | X = [] |
---|
1330 | okey = '' |
---|
1331 | for key in T: |
---|
1332 | if key != okey: X.append(D[key]) #remove duplicate d-spacings |
---|
1333 | okey = key |
---|
1334 | return X |
---|
1335 | |
---|
1336 | def SwapIndx(Axis,H): |
---|
1337 | 'needs doc string' |
---|
1338 | if Axis in [1,-1]: |
---|
1339 | return H |
---|
1340 | elif Axis in [2,-3]: |
---|
1341 | return [H[1],H[2],H[0]] |
---|
1342 | else: |
---|
1343 | return [H[2],H[0],H[1]] |
---|
1344 | |
---|
1345 | def SwapItems(Alist,pos1,pos2): |
---|
1346 | 'exchange 2 items in a list' |
---|
1347 | try: |
---|
1348 | get = Alist[pos1],Alist[pos2] |
---|
1349 | Alist[pos2],Alist[pos1] = get |
---|
1350 | except IndexError: |
---|
1351 | pass |
---|
1352 | return Alist |
---|
1353 | |
---|
1354 | def Rh2Hx(Rh): |
---|
1355 | 'needs doc string' |
---|
1356 | Hx = [0,0,0] |
---|
1357 | Hx[0] = Rh[0]-Rh[1] |
---|
1358 | Hx[1] = Rh[1]-Rh[2] |
---|
1359 | Hx[2] = np.sum(Rh) |
---|
1360 | return Hx |
---|
1361 | |
---|
1362 | def Hx2Rh(Hx): |
---|
1363 | 'needs doc string' |
---|
1364 | Rh = [0,0,0] |
---|
1365 | itk = -Hx[0]+Hx[1]+Hx[2] |
---|
1366 | if itk%3 != 0: |
---|
1367 | return 0 #error - not rhombohedral reflection |
---|
1368 | else: |
---|
1369 | Rh[1] = itk//3 |
---|
1370 | Rh[0] = Rh[1]+Hx[0] |
---|
1371 | Rh[2] = Rh[1]-Hx[1] |
---|
1372 | if Rh[0] < 0: |
---|
1373 | for i in range(3): |
---|
1374 | Rh[i] = -Rh[i] |
---|
1375 | return Rh |
---|
1376 | |
---|
1377 | def CentCheck(Cent,H): |
---|
1378 | 'needs doc string' |
---|
1379 | h,k,l = H |
---|
1380 | if Cent == 'A' and (k+l)%2: |
---|
1381 | return False |
---|
1382 | elif Cent == 'B' and (h+l)%2: |
---|
1383 | return False |
---|
1384 | elif Cent == 'C' and (h+k)%2: |
---|
1385 | return False |
---|
1386 | elif Cent == 'I' and (h+k+l)%2: |
---|
1387 | return False |
---|
1388 | elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2): |
---|
1389 | return False |
---|
1390 | elif Cent == 'R' and (-h+k+l)%3: |
---|
1391 | return False |
---|
1392 | else: |
---|
1393 | return True |
---|
1394 | |
---|
1395 | def RBsymCheck(Atoms,ct,cx,cs,AtLookUp,Amat,RBObjIds,SGData): |
---|
1396 | """ Checks members of a rigid body to see if one is a symmetry equivalent of another. |
---|
1397 | If so the atom site frac is set to zero. |
---|
1398 | |
---|
1399 | :param Atoms: atom array as defined in GSAS-II; modified here |
---|
1400 | :param ct: int location of atom type in Atoms item |
---|
1401 | :param cx: int location of x,y,z,frac in Atoms item |
---|
1402 | :param dict AtLookUp: atom lookup by Id table |
---|
1403 | :param np.array Amat: crystal-to-Cartesian transformation matrix |
---|
1404 | :param list RBObjIds: atom Id belonging to rigid body being tested |
---|
1405 | :param dict SGData: GSAS-II space group info. |
---|
1406 | :returns: Atoms with modified atom frac entries |
---|
1407 | |
---|
1408 | """ |
---|
1409 | for i,Id in enumerate(RBObjIds): |
---|
1410 | XYZo = np.array(Atoms[AtLookUp[Id]][cx:cx+3])%1. |
---|
1411 | typo = Atoms[AtLookUp[Id]][ct] |
---|
1412 | for Jd in RBObjIds[i+1:]: |
---|
1413 | if Atoms[AtLookUp[Jd]][ct] == typo: |
---|
1414 | XYZt = Atoms[AtLookUp[Jd]][cx:cx+3] |
---|
1415 | Xeqv = list(G2spc.GenAtom(np.array(XYZt)%1.,SGData,True)) |
---|
1416 | close = [np.allclose(np.inner(Amat,XYZo),np.inner(Amat,eqv[0]),atol=0.01) for eqv in Xeqv] |
---|
1417 | if True in close: |
---|
1418 | Atoms[AtLookUp[Jd]][cx+3] = 0.0 |
---|
1419 | Sytsym,Mult = G2spc.SytSym(Atoms[AtLookUp[Id]][cx:cx+3],SGData)[:2] |
---|
1420 | Atoms[AtLookUp[Id]][cs] = Sytsym |
---|
1421 | Atoms[AtLookUp[Id]][cs+1] = Mult |
---|
1422 | return Atoms |
---|
1423 | |
---|
1424 | def GetBraviasNum(center,system): |
---|
1425 | """Determine the Bravais lattice number, as used in GenHBravais |
---|
1426 | |
---|
1427 | :param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup) |
---|
1428 | :param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R) |
---|
1429 | 'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup) |
---|
1430 | :return: a number between 0 and 13 |
---|
1431 | or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard) |
---|
1432 | |
---|
1433 | """ |
---|
1434 | if center.upper() == 'F' and system.lower() == 'cubic': |
---|
1435 | return 0 |
---|
1436 | elif center.upper() == 'I' and system.lower() == 'cubic': |
---|
1437 | return 1 |
---|
1438 | elif center.upper() == 'P' and system.lower() == 'cubic': |
---|
1439 | return 2 |
---|
1440 | elif center.upper() == 'R' and system.lower() == 'trigonal': |
---|
1441 | return 3 |
---|
1442 | elif center.upper() == 'P' and system.lower() == 'hexagonal': |
---|
1443 | return 4 |
---|
1444 | elif center.upper() == 'I' and system.lower() == 'tetragonal': |
---|
1445 | return 5 |
---|
1446 | elif center.upper() == 'P' and system.lower() == 'tetragonal': |
---|
1447 | return 6 |
---|
1448 | elif center.upper() == 'F' and system.lower() == 'orthorhombic': |
---|
1449 | return 7 |
---|
1450 | elif center.upper() == 'I' and system.lower() == 'orthorhombic': |
---|
1451 | return 8 |
---|
1452 | elif center.upper() == 'A' and system.lower() == 'orthorhombic': |
---|
1453 | return 9 |
---|
1454 | elif center.upper() == 'B' and system.lower() == 'orthorhombic': |
---|
1455 | return 10 |
---|
1456 | elif center.upper() == 'C' and system.lower() == 'orthorhombic': |
---|
1457 | return 11 |
---|
1458 | elif center.upper() == 'P' and system.lower() == 'orthorhombic': |
---|
1459 | return 12 |
---|
1460 | elif center.upper() == 'C' and system.lower() == 'monoclinic': |
---|
1461 | return 13 |
---|
1462 | elif center.upper() == 'P' and system.lower() == 'monoclinic': |
---|
1463 | return 14 |
---|
1464 | elif center.upper() == 'P' and system.lower() == 'triclinic': |
---|
1465 | return 15 |
---|
1466 | raise ValueError('non-standard Bravais lattice center=%s, cell=%s' % (center,system)) |
---|
1467 | |
---|
1468 | def _GenHBravais_cctbx(dmin, Bravais, A, sg_type, uctbx_unit_cell, miller_index_generator): |
---|
1469 | '''Alternate form of :func:`GenHBravais` that uses CCTBX internals |
---|
1470 | ''' |
---|
1471 | g_inv = np.array([[A[0], A[3]/2, A[4]/2], |
---|
1472 | [A[3]/2, A[1], A[5]/2], |
---|
1473 | [A[4]/2, A[5]/2, A[2]]]) |
---|
1474 | g = np.linalg.inv(g_inv) |
---|
1475 | g_elems = (g[0][0], g[1][1], g[2][2], g[0][1], g[0][2], g[1][2]) |
---|
1476 | try: |
---|
1477 | uc = uctbx_unit_cell(metrical_matrix=g_elems) |
---|
1478 | except ValueError: # this function sometimes receives an A matrix that gives |
---|
1479 | # numbers <0 in the diagonal elems of g. Not sure why. |
---|
1480 | return [] |
---|
1481 | #if sg_type is None: |
---|
1482 | # sg_type = make_sgtype(Bravais) |
---|
1483 | mig = miller_index_generator(uc, sg_type, 0, dmin) |
---|
1484 | result = [] |
---|
1485 | for h,k,l in mig: |
---|
1486 | d = uc.d((h,k,l)) |
---|
1487 | result.append([h, k, l, d, -1]) |
---|
1488 | result.sort(key=lambda l: l[3], reverse=True) |
---|
1489 | return result |
---|
1490 | |
---|
1491 | def GenHBravais(dmin, Bravais, A, cctbx_args=None): |
---|
1492 | """Generate the positionally unique powder diffraction reflections |
---|
1493 | |
---|
1494 | :param dmin: minimum d-spacing in A |
---|
1495 | :param Bravais: lattice type (see GetBraviasNum). Bravais is one of: |
---|
1496 | |
---|
1497 | * 0 F cubic |
---|
1498 | * 1 I cubic |
---|
1499 | * 2 P cubic |
---|
1500 | * 3 R hexagonal (trigonal not rhombohedral) |
---|
1501 | * 4 P hexagonal |
---|
1502 | * 5 I tetragonal |
---|
1503 | * 6 P tetragonal |
---|
1504 | * 7 F orthorhombic |
---|
1505 | * 8 I orthorhombic |
---|
1506 | * 9 A orthorhombic |
---|
1507 | * 10 B orthorhombic |
---|
1508 | * 11 C orthorhombic |
---|
1509 | * 12 P orthorhombic |
---|
1510 | * 13 I monoclinic |
---|
1511 | * 14 A monoclinic |
---|
1512 | * 15 C monoclinic |
---|
1513 | * 16 P monoclinic |
---|
1514 | * 17 P triclinic |
---|
1515 | |
---|
1516 | :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23] |
---|
1517 | :param dict cctbx_args: items defined in CCTBX: |
---|
1518 | |
---|
1519 | * 'sg_type': value from cctbx.sgtbx.space_group_type(symmorphic_sgs[ibrav]) |
---|
1520 | * 'uctbx_unit_cell': pointer to :meth:`cctbx.uctbx.unit_cell` |
---|
1521 | * 'miller_index_generator': pointer to :meth:`cctbx.miller.index_generator` |
---|
1522 | |
---|
1523 | :returns: HKL unique d list of [h,k,l,d,-1] sorted with largest d first |
---|
1524 | |
---|
1525 | """ |
---|
1526 | if cctbx_args: |
---|
1527 | return _GenHBravais_cctbx(dmin, Bravais, A, |
---|
1528 | cctbx_args['sg_type'], cctbx_args['uctbx_unit_cell'], cctbx_args['miller_index_generator']) |
---|
1529 | |
---|
1530 | if Bravais in [9,14]: |
---|
1531 | Cent = 'A' |
---|
1532 | elif Bravais in [10,]: |
---|
1533 | Cent = 'B' |
---|
1534 | elif Bravais in [11,15]: |
---|
1535 | Cent = 'C' |
---|
1536 | elif Bravais in [1,5,8,13]: |
---|
1537 | Cent = 'I' |
---|
1538 | elif Bravais in [0,7]: |
---|
1539 | Cent = 'F' |
---|
1540 | elif Bravais in [3]: |
---|
1541 | Cent = 'R' |
---|
1542 | else: |
---|
1543 | Cent = 'P' |
---|
1544 | Hmax = MaxIndex(dmin,A) |
---|
1545 | dminsq = 1./(dmin**2) |
---|
1546 | HKL = [] |
---|
1547 | if Bravais == 17: #triclinic |
---|
1548 | for l in range(-Hmax[2],Hmax[2]+1): |
---|
1549 | for k in range(-Hmax[1],Hmax[1]+1): |
---|
1550 | hmin = 0 |
---|
1551 | if (k < 0): hmin = 1 |
---|
1552 | if (k ==0 and l < 0): hmin = 1 |
---|
1553 | for h in range(hmin,Hmax[0]+1): |
---|
1554 | H=[h,k,l] |
---|
1555 | rdsq = calc_rDsq(H,A) |
---|
1556 | if 0 < rdsq <= dminsq: |
---|
1557 | HKL.append([h,k,l,rdsq2d(rdsq,6),-1]) |
---|
1558 | elif Bravais in [13,14,15,16]: #monoclinic - b unique |
---|
1559 | Hmax = SwapIndx(2,Hmax) |
---|
1560 | for h in range(Hmax[0]+1): |
---|
1561 | for k in range(-Hmax[1],Hmax[1]+1): |
---|
1562 | lmin = 0 |
---|
1563 | if k < 0:lmin = 1 |
---|
1564 | for l in range(lmin,Hmax[2]+1): |
---|
1565 | [h,k,l] = SwapIndx(-2,[h,k,l]) |
---|
1566 | H = [] |
---|
1567 | if CentCheck(Cent,[h,k,l]): H=[h,k,l] |
---|
1568 | if H: |
---|
1569 | rdsq = calc_rDsq(H,A) |
---|
1570 | if 0 < rdsq <= dminsq: |
---|
1571 | HKL.append([h,k,l,rdsq2d(rdsq,6),-1]) |
---|
1572 | [h,k,l] = SwapIndx(2,[h,k,l]) |
---|
1573 | elif Bravais in [7,8,9,10,11,12]: #orthorhombic |
---|
1574 | for h in range(Hmax[0]+1): |
---|
1575 | for k in range(Hmax[1]+1): |
---|
1576 | for l in range(Hmax[2]+1): |
---|
1577 | H = [] |
---|
1578 | if CentCheck(Cent,[h,k,l]): H=[h,k,l] |
---|
1579 | if H: |
---|
1580 | rdsq = calc_rDsq(H,A) |
---|
1581 | if 0 < rdsq <= dminsq: |
---|
1582 | HKL.append([h,k,l,rdsq2d(rdsq,6),-1]) |
---|
1583 | elif Bravais in [5,6]: #tetragonal |
---|
1584 | for l in range(Hmax[2]+1): |
---|
1585 | for k in range(Hmax[1]+1): |
---|
1586 | for h in range(k,Hmax[0]+1): |
---|
1587 | H = [] |
---|
1588 | if CentCheck(Cent,[h,k,l]): H=[h,k,l] |
---|
1589 | if H: |
---|
1590 | rdsq = calc_rDsq(H,A) |
---|
1591 | if 0 < rdsq <= dminsq: |
---|
1592 | HKL.append([h,k,l,rdsq2d(rdsq,6),-1]) |
---|
1593 | elif Bravais in [3,4]: |
---|
1594 | lmin = 0 |
---|
1595 | if Bravais == 3: lmin = -Hmax[2] #hexagonal/trigonal |
---|
1596 | for l in range(lmin,Hmax[2]+1): |
---|
1597 | for k in range(Hmax[1]+1): |
---|
1598 | hmin = k |
---|
1599 | if l < 0: hmin += 1 |
---|
1600 | for h in range(hmin,Hmax[0]+1): |
---|
1601 | H = [] |
---|
1602 | if CentCheck(Cent,[h,k,l]): H=[h,k,l] |
---|
1603 | if H: |
---|
1604 | rdsq = calc_rDsq(H,A) |
---|
1605 | if 0 < rdsq <= dminsq: |
---|
1606 | HKL.append([h,k,l,rdsq2d(rdsq,6),-1]) |
---|
1607 | |
---|
1608 | else: #cubic |
---|
1609 | for l in range(Hmax[2]+1): |
---|
1610 | for k in range(l,Hmax[1]+1): |
---|
1611 | for h in range(k,Hmax[0]+1): |
---|
1612 | H = [] |
---|
1613 | if CentCheck(Cent,[h,k,l]): H=[h,k,l] |
---|
1614 | if H: |
---|
1615 | rdsq = calc_rDsq(H,A) |
---|
1616 | if 0 < rdsq <= dminsq: |
---|
1617 | HKL.append([h,k,l,rdsq2d(rdsq,6),-1]) |
---|
1618 | return sortHKLd(HKL,True,False) |
---|
1619 | |
---|
1620 | def getHKLmax(dmin,SGData,A): |
---|
1621 | 'finds maximum allowed hkl for given A within dmin' |
---|
1622 | SGLaue = SGData['SGLaue'] |
---|
1623 | if SGLaue in ['3R','3mR']: #Rhombohedral axes |
---|
1624 | Hmax = [0,0,0] |
---|
1625 | cell = A2cell(A) |
---|
1626 | aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3]))) |
---|
1627 | cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3]))) |
---|
1628 | Hmax[0] = Hmax[1] = int(round(aHx/dmin)) |
---|
1629 | Hmax[2] = int(round(cHx/dmin)) |
---|
1630 | #print Hmax,aHx,cHx |
---|
1631 | else: # all others |
---|
1632 | Hmax = MaxIndex(dmin,A) |
---|
1633 | return Hmax |
---|
1634 | |
---|
1635 | def GenHLaue(dmin,SGData,A): |
---|
1636 | """Generate the crystallographically unique powder diffraction reflections |
---|
1637 | for a lattice and Bravais type |
---|
1638 | |
---|
1639 | :param dmin: minimum d-spacing |
---|
1640 | :param SGData: space group dictionary with at least |
---|
1641 | |
---|
1642 | * 'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m' |
---|
1643 | * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F' |
---|
1644 | * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string |
---|
1645 | |
---|
1646 | :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23] |
---|
1647 | :return: HKL = list of [h,k,l,d] sorted with largest d first and is unique |
---|
1648 | part of reciprocal space ignoring anomalous dispersion |
---|
1649 | |
---|
1650 | """ |
---|
1651 | import math |
---|
1652 | SGLaue = SGData['SGLaue'] |
---|
1653 | SGLatt = SGData['SGLatt'] |
---|
1654 | SGUniq = SGData['SGUniq'] |
---|
1655 | #finds maximum allowed hkl for given A within dmin |
---|
1656 | Hmax = getHKLmax(dmin,SGData,A) |
---|
1657 | |
---|
1658 | dminsq = 1./(dmin**2) |
---|
1659 | HKL = [] |
---|
1660 | if SGLaue == '-1': #triclinic |
---|
1661 | for l in range(-Hmax[2],Hmax[2]+1): |
---|
1662 | for k in range(-Hmax[1],Hmax[1]+1): |
---|
1663 | hmin = 0 |
---|
1664 | if (k < 0) or (k ==0 and l < 0): hmin = 1 |
---|
1665 | for h in range(hmin,Hmax[0]+1): |
---|
1666 | H = [] |
---|
1667 | if CentCheck(SGLatt,[h,k,l]): H=[h,k,l] |
---|
1668 | if H: |
---|
1669 | rdsq = calc_rDsq(H,A) |
---|
1670 | if 0 < rdsq <= dminsq: |
---|
1671 | HKL.append([h,k,l,1./math.sqrt(rdsq)]) |
---|
1672 | elif SGLaue == '2/m': #monoclinic |
---|
1673 | axisnum = 1 + ['a','b','c'].index(SGUniq) |
---|
1674 | Hmax = SwapIndx(axisnum,Hmax) |
---|
1675 | for h in range(Hmax[0]+1): |
---|
1676 | for k in range(-Hmax[1],Hmax[1]+1): |
---|
1677 | lmin = 0 |
---|
1678 | if k < 0:lmin = 1 |
---|
1679 | for l in range(lmin,Hmax[2]+1): |
---|
1680 | [h,k,l] = SwapIndx(-axisnum,[h,k,l]) |
---|
1681 | H = [] |
---|
1682 | if CentCheck(SGLatt,[h,k,l]): H=[h,k,l] |
---|
1683 | if H: |
---|
1684 | rdsq = calc_rDsq(H,A) |
---|
1685 | if 0 < rdsq <= dminsq: |
---|
1686 | HKL.append([h,k,l,1./math.sqrt(rdsq)]) |
---|
1687 | [h,k,l] = SwapIndx(axisnum,[h,k,l]) |
---|
1688 | elif SGLaue in ['mmm','4/m','6/m']: #orthorhombic |
---|
1689 | for l in range(Hmax[2]+1): |
---|
1690 | for h in range(Hmax[0]+1): |
---|
1691 | kmin = 1 |
---|
1692 | if SGLaue == 'mmm' or h ==0: kmin = 0 |
---|
1693 | for k in range(kmin,Hmax[1]+1): |
---|
1694 | H = [] |
---|
1695 | if CentCheck(SGLatt,[h,k,l]): H=[h,k,l] |
---|
1696 | if H: |
---|
1697 | rdsq = calc_rDsq(H,A) |
---|
1698 | if 0 < rdsq <= dminsq: |
---|
1699 | HKL.append([h,k,l,1./math.sqrt(rdsq)]) |
---|
1700 | elif SGLaue in ['4/mmm','6/mmm']: #tetragonal & hexagonal |
---|
1701 | for l in range(Hmax[2]+1): |
---|
1702 | for h in range(Hmax[0]+1): |
---|
1703 | for k in range(h+1): |
---|
1704 | H = [] |
---|
1705 | if CentCheck(SGLatt,[h,k,l]): H=[h,k,l] |
---|
1706 | if H: |
---|
1707 | rdsq = calc_rDsq(H,A) |
---|
1708 | if 0 < rdsq <= dminsq: |
---|
1709 | HKL.append([h,k,l,1./math.sqrt(rdsq)]) |
---|
1710 | elif SGLaue in ['3m1','31m','3','3R','3mR']: #trigonals |
---|
1711 | for l in range(-Hmax[2],Hmax[2]+1): |
---|
1712 | hmin = 0 |
---|
1713 | if l < 0: hmin = 1 |
---|
1714 | for h in range(hmin,Hmax[0]+1): |
---|
1715 | if SGLaue in ['3R','3']: |
---|
1716 | kmax = h |
---|
1717 | kmin = -int((h-1.)/2.) |
---|
1718 | else: |
---|
1719 | kmin = 0 |
---|
1720 | kmax = h |
---|
1721 | if SGLaue in ['3m1','3mR'] and l < 0: kmax = h-1 |
---|
1722 | if SGLaue == '31m' and l < 0: kmin = 1 |
---|
1723 | for k in range(kmin,kmax+1): |
---|
1724 | H = [] |
---|
1725 | if CentCheck(SGLatt,[h,k,l]): H=[h,k,l] |
---|
1726 | if SGLaue in ['3R','3mR']: |
---|
1727 | H = Hx2Rh(H) |
---|
1728 | if H: |
---|
1729 | rdsq = calc_rDsq(H,A) |
---|
1730 | if 0 < rdsq <= dminsq: |
---|
1731 | HKL.append([H[0],H[1],H[2],1./math.sqrt(rdsq)]) |
---|
1732 | else: #cubic |
---|
1733 | for h in range(Hmax[0]+1): |
---|
1734 | for k in range(h+1): |
---|
1735 | lmin = 0 |
---|
1736 | lmax = k |
---|
1737 | if SGLaue =='m3': |
---|
1738 | lmax = h-1 |
---|
1739 | if h == k: lmax += 1 |
---|
1740 | for l in range(lmin,lmax+1): |
---|
1741 | H = [] |
---|
1742 | if CentCheck(SGLatt,[h,k,l]): H=[h,k,l] |
---|
1743 | if H: |
---|
1744 | rdsq = calc_rDsq(H,A) |
---|
1745 | if 0 < rdsq <= dminsq: |
---|
1746 | HKL.append([h,k,l,1./math.sqrt(rdsq)]) |
---|
1747 | return sortHKLd(HKL,True,True) |
---|
1748 | |
---|
1749 | def GenPfHKLs(nMax,SGData,A): |
---|
1750 | """Generate the unique pole figure reflections for a lattice and Bravais type. |
---|
1751 | Min d-spacing=1.0A & no more than nMax returned |
---|
1752 | |
---|
1753 | :param nMax: maximum number of hkls returned |
---|
1754 | :param SGData: space group dictionary with at least |
---|
1755 | |
---|
1756 | * 'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m' |
---|
1757 | * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F' |
---|
1758 | * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string |
---|
1759 | |
---|
1760 | :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23] |
---|
1761 | :return: HKL = list of 'h k l' strings sorted with largest d first; no duplicate zones |
---|
1762 | |
---|
1763 | """ |
---|
1764 | HKL = np.array(GenHLaue(1.0,SGData,A)).T[:3].T #strip d-spacings |
---|
1765 | N = min(nMax,len(HKL)) |
---|
1766 | return ['%d %d %d'%(h[0],h[1],h[2]) for h in HKL[:N]] |
---|
1767 | |
---|
1768 | def GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A): |
---|
1769 | 'needs a doc string' |
---|
1770 | ifMag = False |
---|
1771 | if 'MagSpGrp' in SGData: |
---|
1772 | ifMag = True |
---|
1773 | HKLs = [] |
---|
1774 | vec = np.array(Vec) |
---|
1775 | vstar = np.sqrt(calc_rDsq(vec,A)) #find extra needed for -n SS reflections |
---|
1776 | dvec = 1./(maxH*vstar+1./dmin) |
---|
1777 | HKL = GenHLaue(dvec,SGData,A) |
---|
1778 | SSdH = [vec*h for h in range(-maxH,maxH+1)] |
---|
1779 | SSdH = dict(zip(range(-maxH,maxH+1),SSdH)) |
---|
1780 | for h,k,l,d in HKL: |
---|
1781 | ext = G2spc.GenHKLf([h,k,l],SGData)[0] #h,k,l must be integral values here |
---|
1782 | if not ext and d >= dmin: |
---|
1783 | HKLs.append([h,k,l,0,d]) |
---|
1784 | for dH in SSdH: |
---|
1785 | if dH: |
---|
1786 | DH = SSdH[dH] |
---|
1787 | H = [h+DH[0],k+DH[1],l+DH[2]] |
---|
1788 | d = 1./np.sqrt(calc_rDsq(H,A)) |
---|
1789 | if d >= dmin: |
---|
1790 | HKLM = np.array([h,k,l,dH]) |
---|
1791 | if (G2spc.checkSSLaue([h,k,l,dH],SGData,SSGData) and G2spc.checkSSextc(HKLM,SSGData)) or ifMag: |
---|
1792 | HKLs.append([h,k,l,dH,d]) |
---|
1793 | return HKLs |
---|
1794 | |
---|
1795 | def LaueUnique2(SGData,refList): |
---|
1796 | ''' Impose Laue symmetry on hkl |
---|
1797 | |
---|
1798 | :param SGData: space group data from 'P '+Laue |
---|
1799 | :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted |
---|
1800 | |
---|
1801 | :return: HKLF new reflection array with imposed Laue symmetry |
---|
1802 | ''' |
---|
1803 | for ref in refList: |
---|
1804 | H = ref[:3] |
---|
1805 | Uniq = G2spc.GenHKLf(H,SGData)[2] |
---|
1806 | Uniq = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Uniq,2),1),0) |
---|
1807 | ref[:3] = Uniq[-1] |
---|
1808 | return refList |
---|
1809 | |
---|
1810 | def LaueUnique(Laue,HKLF): |
---|
1811 | ''' Impose Laue symmetry on hkl |
---|
1812 | |
---|
1813 | :param str Laue: Laue symbol, as below |
---|
1814 | |
---|
1815 | centrosymmetric Laue groups:: |
---|
1816 | |
---|
1817 | ['-1','2/m','112/m','2/m11','mmm','-42m','-4m2','4/mmm','-3','-3m', |
---|
1818 | '-31m','-3m1','6/m','6/mmm','m3','m3m'] |
---|
1819 | |
---|
1820 | noncentrosymmetric Laue groups:: |
---|
1821 | |
---|
1822 | ['1','2','211','112','m','m11','11m','222','mm2','m2m','2mm', |
---|
1823 | '4','-4','422','4mm','3','312','321','3m','31m','3m1','6','-6', |
---|
1824 | '622','6mm','-62m','-6m2','23','432','-43m'] |
---|
1825 | |
---|
1826 | :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted |
---|
1827 | |
---|
1828 | :returns: HKLF new reflection array with imposed Laue symmetry |
---|
1829 | ''' |
---|
1830 | |
---|
1831 | HKLFT = HKLF.T |
---|
1832 | mat41 = np.array([[0,1,0],[-1,0,0],[0,0,1]]) #hkl -> k,-h,l |
---|
1833 | mat43 = np.array([[0,-1,0],[1,0,0],[0,0,1]]) #hkl -> -k,h,l |
---|
1834 | mat4bar = np.array([[0,-1,0],[1,0,0],[0,0,-1]]) #hkl -> k,-h,-l |
---|
1835 | mat31 = np.array([[-1,-1,0],[1,0,0],[0,0,1]]) #hkl -> ihl = -h-k,h,l |
---|
1836 | mat32 = np.array([[0,1,0],[-1,-1,0],[0,0,1]]) #hkl -> kil = k,-h-k,l |
---|
1837 | matd3 = np.array([[0,1,0],[0,0,1],[1,0,0]]) #hkl -> k,l,h |
---|
1838 | matd3q = np.array([[0,0,-1],[-1,0,0],[0,1,0]]) #hkl -> -l,-h,k |
---|
1839 | matd3t = np.array([[0,0,-1],[1,0,0],[0,-1,0]]) #hkl -> -l,h,-k |
---|
1840 | mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]]) #hkl -> h+k,-h,l really 65 |
---|
1841 | matdm = np.array([[0,1,0],[1,0,0],[0,0,1]]) #hkl -> k,h,l |
---|
1842 | matdmp = np.array([[-1,-1,0],[0,1,0],[0,0,1]]) #hkl -> -h-k,k,l |
---|
1843 | matkm = np.array([[-1,0,0],[1,1,0],[0,0,1]]) #hkl -> -h,h+k,l |
---|
1844 | matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]]) #hkl -> k,h,-l |
---|
1845 | matdm3 = np.array([[1,0,0],[0,0,1],[0,1,0]]) #hkl -> h,l,k |
---|
1846 | mat2d43 = np.array([[0,1,0],[1,0,0],[0,0,1]]) #hkl -> k,-h,l |
---|
1847 | matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]]) #hkl -> -h,-i,-l |
---|
1848 | #triclinic |
---|
1849 | if Laue == '1': #ok |
---|
1850 | pass |
---|
1851 | elif Laue == '-1': #ok |
---|
1852 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
1853 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
1854 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
1855 | #monoclinic |
---|
1856 | #noncentrosymmetric - all ok |
---|
1857 | elif Laue == '2': |
---|
1858 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1859 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1860 | elif Laue == '1 1 2': |
---|
1861 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1862 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1863 | elif Laue == '2 1 1': |
---|
1864 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
1865 | HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
1866 | elif Laue == 'm': |
---|
1867 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1868 | elif Laue == 'm 1 1': |
---|
1869 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3]) |
---|
1870 | elif Laue == '1 1 m': |
---|
1871 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1872 | #centrosymmetric - all ok |
---|
1873 | elif Laue == '2/m 1 1': |
---|
1874 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
1875 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3]) |
---|
1876 | HKLFT[:3] = np.where((HKLFT[2]*HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1877 | elif Laue == '2/m': |
---|
1878 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1879 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1880 | HKLFT[:3] = np.where((HKLFT[0]*HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1881 | elif Laue == '1 1 2/m': |
---|
1882 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1883 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1884 | HKLFT[:3] = np.where((HKLFT[1]*HKLFT[2]==0)&(HKLFT[0]<0),HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3]) |
---|
1885 | #orthorhombic |
---|
1886 | #noncentrosymmetric - all OK |
---|
1887 | elif Laue == '2 2 2': |
---|
1888 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1889 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
1890 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1891 | HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1892 | elif Laue == 'm m 2': |
---|
1893 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3]) |
---|
1894 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1895 | elif Laue == '2 m m': |
---|
1896 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1897 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1898 | elif Laue == 'm 2 m': |
---|
1899 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3]) |
---|
1900 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1901 | #centrosymmetric - all ok |
---|
1902 | elif Laue == 'm m m': |
---|
1903 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3]) |
---|
1904 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1905 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1906 | #tetragonal |
---|
1907 | #noncentrosymmetric - all ok |
---|
1908 | elif Laue == '4': |
---|
1909 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1910 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3]) |
---|
1911 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3]) |
---|
1912 | elif Laue == '-4': |
---|
1913 | HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1914 | HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3]) |
---|
1915 | HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1916 | HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3]) |
---|
1917 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1918 | elif Laue == '4 2 2': |
---|
1919 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
1920 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1921 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3]) |
---|
1922 | HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3]) |
---|
1923 | HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) #in lieu od 2-fold |
---|
1924 | elif Laue == '4 m m': |
---|
1925 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1926 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1927 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3]) |
---|
1928 | HKLFT[:3] = np.where(HKLFT[0]<HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) |
---|
1929 | elif Laue == '-4 2 m': |
---|
1930 | HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1931 | HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3]) |
---|
1932 | HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1933 | HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3]) |
---|
1934 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1935 | HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) |
---|
1936 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1937 | elif Laue == '-4 m 2': |
---|
1938 | HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3]) |
---|
1939 | HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1940 | HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<=0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3]) |
---|
1941 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1942 | HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]==0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3]) |
---|
1943 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1944 | HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) |
---|
1945 | #centrosymmetric - all ok |
---|
1946 | elif Laue == '4/m': |
---|
1947 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1948 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1949 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3]) |
---|
1950 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3]) |
---|
1951 | elif Laue == '4/m m m': |
---|
1952 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1953 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1954 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3]) |
---|
1955 | HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3]) |
---|
1956 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) |
---|
1957 | #trigonal - all hex cell |
---|
1958 | #noncentrosymmetric - all ok |
---|
1959 | elif Laue == '3': |
---|
1960 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1961 | HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1962 | HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3]) |
---|
1963 | elif Laue == '3 1 2': |
---|
1964 | HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3]) |
---|
1965 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1966 | HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1967 | HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3]) |
---|
1968 | HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3]) |
---|
1969 | elif Laue == '3 2 1': |
---|
1970 | HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1971 | HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1972 | HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3]) |
---|
1973 | HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3]) |
---|
1974 | HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T |
---|
1975 | HKLFT[:3] = np.where((HKLFT[0]!=0)&(HKLFT[2]>0)&(HKLFT[0]==-2*HKLFT[1]),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
1976 | elif Laue == '3 1 m': |
---|
1977 | HKLFT[:3] = np.where(HKLFT[0]>=HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1978 | HKLFT[:3] = np.where(2*HKLFT[1]<-HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1979 | HKLFT[:3] = np.where(HKLFT[1]>-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdmp[nxs,:,:])).T,HKLFT[:3]) |
---|
1980 | HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T |
---|
1981 | elif (Laue == '3 m 1' or Laue == '3 m'): |
---|
1982 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1983 | HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1984 | HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3]) |
---|
1985 | #centrosymmetric |
---|
1986 | elif Laue == '-3': #ok |
---|
1987 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
1988 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1989 | HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1990 | HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3]) |
---|
1991 | HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]<0),-np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3]) |
---|
1992 | HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3]) |
---|
1993 | elif (Laue == '-3 m 1' or Laue == '-3 m'): #ok |
---|
1994 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1995 | HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
1996 | HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3]) |
---|
1997 | HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3]) |
---|
1998 | HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3]) |
---|
1999 | elif Laue == '-3 1 m': #ok |
---|
2000 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
2001 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2002 | HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2003 | HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3]) |
---|
2004 | HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3]) |
---|
2005 | HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) |
---|
2006 | #hexagonal |
---|
2007 | #noncentrosymmetric |
---|
2008 | elif Laue == '6': #ok |
---|
2009 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2010 | HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2011 | HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3]) |
---|
2012 | HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3]) |
---|
2013 | elif Laue == '-6': #ok |
---|
2014 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2015 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2016 | HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2017 | HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3]) |
---|
2018 | elif Laue == '6 2 2': #ok |
---|
2019 | HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3]) |
---|
2020 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2021 | HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2022 | HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3]) |
---|
2023 | HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) |
---|
2024 | HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) |
---|
2025 | elif Laue == '6 m m': #ok |
---|
2026 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2027 | HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2028 | HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3]) |
---|
2029 | HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3]) |
---|
2030 | HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) |
---|
2031 | elif Laue == '-6 m 2': #ok |
---|
2032 | HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3]) |
---|
2033 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2034 | HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2035 | HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3]) |
---|
2036 | HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3]) |
---|
2037 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2038 | elif Laue == '-6 2 m': #ok |
---|
2039 | HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3]) |
---|
2040 | HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2041 | HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2042 | HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3]) |
---|
2043 | HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3]) |
---|
2044 | HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T |
---|
2045 | HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3]) |
---|
2046 | HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) |
---|
2047 | #centrosymmetric |
---|
2048 | elif Laue == '6/m': #ok |
---|
2049 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2050 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2051 | HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2052 | HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3]) |
---|
2053 | HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3]) |
---|
2054 | elif Laue == '6/m m m': #ok |
---|
2055 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2056 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2057 | HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3]) |
---|
2058 | HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3]) |
---|
2059 | HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm.T[nxs,:,:])).T,HKLFT[:3]) |
---|
2060 | #cubic - all ok |
---|
2061 | #noncentrosymmetric - |
---|
2062 | elif Laue == '2 3': |
---|
2063 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
2064 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
2065 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2066 | HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2067 | HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3]) |
---|
2068 | HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3]) |
---|
2069 | HKLFT[:3] = np.where((HKLFT[2]<0)&((HKLFT[0]>-HKLFT[2])|(HKLFT[1]>-HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3t[nxs,:,:])).T,HKLFT[:3]) |
---|
2070 | HKLFT[:3] = np.where((HKLFT[2]<0)&((HKLFT[0]>-HKLFT[2])|(HKLFT[1]>=-HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3t[nxs,:,:])).T,HKLFT[:3]) |
---|
2071 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2072 | elif Laue == '4 3 2': |
---|
2073 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3]) |
---|
2074 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
2075 | HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3]) |
---|
2076 | HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3]) |
---|
2077 | HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) #in lieu od 2-fold |
---|
2078 | HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3]) |
---|
2079 | HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3]) |
---|
2080 | HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat2d43[nxs,:,:])).T,HKLFT[:3]) |
---|
2081 | elif Laue == '-4 3 m': |
---|
2082 | HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
2083 | HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3]) |
---|
2084 | HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3]) |
---|
2085 | HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3]) |
---|
2086 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2087 | HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) |
---|
2088 | HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2089 | HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3]) |
---|
2090 | HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3]) |
---|
2091 | HKLFT[:3] = np.where((HKLFT[2]>=0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) |
---|
2092 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2093 | HKLFT[:3] = np.where((HKLFT[0]<0)&(HKLFT[2]<-HKLFT[0])&(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3q[nxs,:,:])).T,HKLFT[:3]) |
---|
2094 | HKLFT[:3] = np.where((HKLFT[0]<0)&(HKLFT[2]>=-HKLFT[0])&(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matdm3[nxs,:,:])).T,HKLFT[:3]) |
---|
2095 | #centrosymmetric |
---|
2096 | elif Laue == 'm 3': |
---|
2097 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3]) |
---|
2098 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) |
---|
2099 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2100 | HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3]) |
---|
2101 | HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3]) |
---|
2102 | elif Laue == 'm 3 m': |
---|
2103 | HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3]) |
---|
2104 | HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) |
---|
2105 | HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3]) |
---|
2106 | HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3]) |
---|
2107 | HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3]) |
---|
2108 | HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3]) |
---|
2109 | return HKLFT.T |
---|
2110 | |
---|
2111 | #Spherical harmonics routines |
---|
2112 | def RBChk(sytsym,L,M): |
---|
2113 | '''finds symmetry rules for spherical harmonic coefficients for site symmetries |
---|
2114 | :param str SGLaue: Laue symbol or sytsym symbol |
---|
2115 | :param int L: principal harmonic term |
---|
2116 | :param int M: second harmonic term; can be -L <= M <= L |
---|
2117 | :returns True if allowed and sign for term |
---|
2118 | NB: not complete for all possible site symmetries! Many are missing |
---|
2119 | Based on Tables 2 & 4 of M. Kara & K. Kurki-Suonio, Acta Cryst. A37, 201-210 (1981). |
---|
2120 | ''' |
---|
2121 | if M <= L: |
---|
2122 | if sytsym == '53m': |
---|
2123 | if not L%2 and M > 0: |
---|
2124 | if L in [6,10,12,16,18]: |
---|
2125 | if L%12 == 2: |
---|
2126 | if M <= L//12: return True,1.0 |
---|
2127 | else: |
---|
2128 | if M <= L//12+1: return True,1.0 |
---|
2129 | elif sytsym == '23': #cubics use different Fourier expansion than those below |
---|
2130 | if 2 < L < 11 and [L,M] in [[3,1],[4,1],[6,1],[6,2],[7,1],[8,1],[9,1],[9,2],[10,1],[10,2]]: |
---|
2131 | return True,1.0 |
---|
2132 | elif sytsym == 'm3': |
---|
2133 | if 2 < L < 11 and [L,M] in [[4,1],[6,1],[6,2],[8,1],[10,1],[10,2]]: |
---|
2134 | return True,1.0 |
---|
2135 | elif sytsym == '432': |
---|
2136 | if 2 < L < 11 and [L,M] in [[4,1],[6,1],[8,1],[9,2],[10,1]]: |
---|
2137 | return True,1.0 |
---|
2138 | elif sytsym == '-43m': |
---|
2139 | if 2 < L < 11 and [L,M] in [[3,1],[4,1],[6,1],[7,1],[8,1],[9,1],[10,1]]: |
---|
2140 | return True,1.0 |
---|
2141 | elif sytsym == 'm3m': #correct for L < 21 by generator |
---|
2142 | if not L%2 and M > 0: |
---|
2143 | if L%12 == 2: |
---|
2144 | if M <= L//12: return True,1.0 |
---|
2145 | else: |
---|
2146 | if M <= L//12+1: return True,1.0 |
---|
2147 | elif sytsym == '6': |
---|
2148 | if not M%6: return True,1.0 #P? |
---|
2149 | elif sytsym == '-6': #L=M+2J |
---|
2150 | if L != 1 and not M%3: #P? |
---|
2151 | if not L%2 and not M%6: return True,1.0 |
---|
2152 | elif L%2 and (M//3)%2: return True,1.0 |
---|
2153 | elif sytsym == '6/m': |
---|
2154 | if not L%2 and not M%6: return True,1.0 #P? |
---|
2155 | elif sytsym == '622': |
---|
2156 | if not M%6: return True,-1.**M |
---|
2157 | elif sytsym == '6mm': |
---|
2158 | if not M%6: return True,1.0 |
---|
2159 | elif sytsym in ['-6m2(100)','-6m2']: #L=M+2J |
---|
2160 | if L != 1 and not M%3: |
---|
2161 | if not L%2 and not M%6: return True,1.0 |
---|
2162 | elif L%2 and (M//3)%2: return True,1.0 |
---|
2163 | elif sytsym == '-6m2(120)': #L=M+2J |
---|
2164 | if L != 1 and not M%3: |
---|
2165 | if not L%2 and not M%6: return True,1.0 |
---|
2166 | elif L%2 and (M//3)%2: return True,-1.**M |
---|
2167 | elif sytsym == '6/mmm': |
---|
2168 | if not L%2 and not M%6: return True,1.0 |
---|
2169 | elif sytsym == '4(z)': |
---|
2170 | if not M%4: return True,1.0 #P? |
---|
2171 | elif sytsym == '-4(z)': #m=2l-4j |
---|
2172 | if L%2 and (M//2)%2: return True,1.0 #P? |
---|
2173 | if not L%2 and not (M//2)%2: return True,1.0 |
---|
2174 | elif sytsym == '4/m(z)': |
---|
2175 | if not M%4: return True,1.0 #P? |
---|
2176 | elif sytsym == '422(z)': |
---|
2177 | if not M%4: return True,-1.0**L |
---|
2178 | elif sytsym == '4mm(z)': |
---|
2179 | if not M%4: return True,1.0 |
---|
2180 | elif sytsym in ['-42m(z)','-42m']: #m=2l-4j |
---|
2181 | if L%2 and (M//2)%2: return True,1.0 |
---|
2182 | if not L%2 and not (M//2)%2: return True,-1.0**L |
---|
2183 | elif sytsym == '-4m2(z)': #m=2l-4j |
---|
2184 | if L%2 and (M//2)%2: return True,1.0 |
---|
2185 | if not L%2 and not (M//2)%2: return True,1.0 |
---|
2186 | elif sytsym == '4/mmm(z)': |
---|
2187 | if not L%2 and not M%4: return True,1.0 |
---|
2188 | elif sytsym == '3' or sytsym == '3(111)': |
---|
2189 | if not M%3: return True,1.0 #P? |
---|
2190 | elif sytsym == '-3' or sytsym == '-3(111)': |
---|
2191 | if not L%2 and not M%3: return True,1.0 #P? |
---|
2192 | elif sytsym == '32(100)' or sytsym == '32(111)': |
---|
2193 | if not M%3: return True,-1.0**L |
---|
2194 | elif sytsym == '32(120)': |
---|
2195 | if not M%3: return True,-1.0**(L-M) |
---|
2196 | elif sytsym in ['3m(100)','3m(111)']: |
---|
2197 | if not M%3: return True,-1.0**M |
---|
2198 | elif sytsym == '3m(120)': |
---|
2199 | if not M%3: return True,1.0 |
---|
2200 | elif sytsym in ['-3m(100)','-3m(111)','-3m']: |
---|
2201 | if not L%2 and not M%3: return True,-1.0**M |
---|
2202 | elif sytsym == '-3m(120)': |
---|
2203 | if not L%2 and not M%3: return True,1.0 |
---|
2204 | elif '222' in sytsym: |
---|
2205 | if M%2: return True,-1.0**L |
---|
2206 | elif 'mm2(x)' in sytsym: #m=l-2j |
---|
2207 | if L%2 and M%2: return True,1.0 #both odd |
---|
2208 | if not L%2 and not M%2: return True,1.0 #both even |
---|
2209 | elif 'mm2(y)' in sytsym: #m=l-2j |
---|
2210 | if L%2 and M%2: return True,-1.0**L #both odd |
---|
2211 | if not L%2 and not M%2: return True,-1.0**L #both even |
---|
2212 | elif 'mm2(z)' in sytsym: |
---|
2213 | if M%2: return True,1.0 |
---|
2214 | elif 'mmm' in sytsym : |
---|
2215 | if not L%2 and not M%2: return True,1.0 |
---|
2216 | elif sytsym == '2(x)': |
---|
2217 | return True,-1.0**(L-M) |
---|
2218 | elif sytsym == '2(y)': |
---|
2219 | return True,-1.0**L |
---|
2220 | elif sytsym == '2(z)': |
---|
2221 | if not M%2: return True,1.0 #P? |
---|
2222 | elif sytsym == 'm(x)': |
---|
2223 | if not L%2 : return True,-1.0**M |
---|
2224 | elif sytsym == 'm(y)': |
---|
2225 | return True,1.0 |
---|
2226 | elif sytsym == 'm(z)': #m=l-2j |
---|
2227 | if L%2 and M%2: return True,1.0 #P? |
---|
2228 | if not L%2 and not M%2: return True,1.0 #P? |
---|
2229 | elif sytsym == '2/m(x)': |
---|
2230 | if not L%2 : return True,-1.0**M |
---|
2231 | elif sytsym in ['2/m(y)','2/m']: |
---|
2232 | if not L%2: return True,1.0 |
---|
2233 | elif sytsym == '2/m(z)': |
---|
2234 | if not L%2 and not M%2: return True,1.0 |
---|
2235 | elif sytsym == '1': #P? |
---|
2236 | return True,1.0 |
---|
2237 | elif sytsym == '-1': #P? |
---|
2238 | if not L%2: return True,1.0 |
---|
2239 | return False,0. |
---|
2240 | |
---|
2241 | def RBsymChk(RBsym,cubic,coefNames,L=18): |
---|
2242 | '''imposes rigid body symmetry on spherical harmonics terms |
---|
2243 | Key problem is noncubic RB symmetries in cubic site symmetries & vice versa. |
---|
2244 | ''' |
---|
2245 | # cubicsigns = {'C(3,1)c':[-1,],'C(4,1)c':[1,1,],'C(6,1)c':[1,-1],'C(6,2)c':[1,-1],'C(7,1)c':[1,-1],'C(8,1)c':[1,1,1], |
---|
2246 | # 'C(9,1)c':[1,-1],'C(9,2)c':[1,-1],'C(10,1)c':[1,-1,-1],'C(10,2)c':[1,1,-1]} |
---|
2247 | # cubicnames = {'C(3,1)c':['C(3,2)',],'C(4,1)c':['C(4,0)','C(4,4)'],'C(6,1)c':['C(6,0)','C(6,4)'], |
---|
2248 | # 'C(6,2)c':['C(6,2)','C(6,6)'],'C(7,1)c':['C(7,2)','C(7,6)'],'C(8,1)c':['C(8,0)','C(8,4)','C(8,8)'], |
---|
2249 | # 'C(9,1)c':['C(9,2)','C((9,6)'],'C(9,2)c':['C(9,4)','C(9,8)'], |
---|
2250 | # 'C(10,1)c':['C(10,0)','C(10,4)','C(10,8)'],'C(10,2)c':['C(10,2)','C(10,6)','C(10,10)']} |
---|
2251 | newNames = [] |
---|
2252 | newSgns = [] |
---|
2253 | if cubic: #sytsym is a cubic site |
---|
2254 | if RBsym in ['53m','532']: |
---|
2255 | for name in coefNames: |
---|
2256 | LM = eval(name[1:-1]) |
---|
2257 | if LM[0] in [6,10,12,16,18]: |
---|
2258 | newNames.append(name) |
---|
2259 | newSgns.append(1.0) |
---|
2260 | elif RBsym in ['m3m','-43m']: #take all terms? |
---|
2261 | for name in coefNames: |
---|
2262 | LM = eval(name[1:-1]) |
---|
2263 | rbChk,sgn = RBChk(RBsym,LM[0],LM[1]) |
---|
2264 | if rbChk: |
---|
2265 | newNames.append(name) |
---|
2266 | newSgns.append(1.0) |
---|
2267 | else: #RBsym not cubic or icosahedral |
---|
2268 | for name in coefNames: #these are cubic names |
---|
2269 | LM = eval(name[1:-1]) |
---|
2270 | if (LM[0]+LM[1])%2: #even L odd M or vv |
---|
2271 | if LM[0]%2: |
---|
2272 | M = [4*m for m in range(LM[0]//2)[1:] if 4*m <= LM[0]] |
---|
2273 | else: |
---|
2274 | M = [4*m for m in range(LM[0]//2) if 4*m <= LM[0]] |
---|
2275 | else: #both even or both odd |
---|
2276 | M = [4*m+2 for m in range(LM[0]//2) if 4*m+2 <= LM[0]] |
---|
2277 | for m in M: |
---|
2278 | rbChk,sgn = RBChk(RBsym,LM[0],m) |
---|
2279 | if rbChk: |
---|
2280 | newNames.append('C(%d,%d)'%(LM[0],m)) |
---|
2281 | newSgns.append(sgn) |
---|
2282 | else: |
---|
2283 | if RBsym in ['m3m','-43m','53m']: #force mol. sym. here |
---|
2284 | for L in range(L+1): |
---|
2285 | cubNames,cubSgns = GenShCoeff(RBsym,L) |
---|
2286 | newNames += cubNames |
---|
2287 | newSgns += cubSgns |
---|
2288 | else: |
---|
2289 | for name in coefNames: |
---|
2290 | LM = eval(name[1:]) |
---|
2291 | rbChk,sgn = RBChk(RBsym,LM[0],LM[1]) |
---|
2292 | if rbChk: |
---|
2293 | newNames.append(name) |
---|
2294 | newSgns.append(sgn) |
---|
2295 | return newNames,newSgns |
---|
2296 | |
---|
2297 | def GenRBCoeff(sytsym,RBsym,L): |
---|
2298 | '''imposes rigid body symmetry on spherical harmonics terms |
---|
2299 | Key problem is noncubic RB symmetries in cubic site symmetries & vice versa. |
---|
2300 | ''' |
---|
2301 | coefNames = [] |
---|
2302 | coefSgns = [] |
---|
2303 | cubic = False |
---|
2304 | if sytsym in ['23','m3','432','-43m','m3m']: |
---|
2305 | cubic = True |
---|
2306 | for iord in range(L+1): |
---|
2307 | if not iord: continue |
---|
2308 | for n in range(iord+1): |
---|
2309 | rbChk,sgn = RBChk(sytsym,iord,n) |
---|
2310 | if rbChk: |
---|
2311 | if cubic: |
---|
2312 | coefNames.append('C(%d,%d)c'%(iord,n)) |
---|
2313 | else: |
---|
2314 | coefNames.append('C(%d,%d)'%(iord,n)) |
---|
2315 | coefSgns.append(sgn) |
---|
2316 | if RBsym == '1': |
---|
2317 | return coefNames,coefSgns |
---|
2318 | newNames,newSgns = RBsymChk(RBsym,cubic,coefNames,L) |
---|
2319 | return newNames,newSgns |
---|
2320 | |
---|
2321 | def GenShCoeff(sytsym,L): |
---|
2322 | coefNames = [] |
---|
2323 | coefSgns = [] |
---|
2324 | cubic = False |
---|
2325 | if sytsym in ['23','m3','432','-43m','m3m','53m']: |
---|
2326 | cubic = True |
---|
2327 | for n in range(L+1): |
---|
2328 | rbChk,sgn = RBChk(sytsym,L,n) |
---|
2329 | if rbChk: |
---|
2330 | if cubic: |
---|
2331 | coefNames.append('C(%d,%d)c'%(L,n)) |
---|
2332 | else: |
---|
2333 | coefNames.append('C(%d,%d)'%(L,n)) |
---|
2334 | coefSgns.append(sgn) |
---|
2335 | newNames,newSgns = RBsymChk(sytsym,cubic,coefNames,L) |
---|
2336 | return newNames,newSgns |
---|
2337 | |
---|
2338 | def OdfChk(SGLaue,L,M): |
---|
2339 | '''finds symmetry rules for spherical harmonic coefficients for Laue groups |
---|
2340 | :param str SGLaue: Laue symbol |
---|
2341 | :param int L: principal harmonic term; only evens are used |
---|
2342 | :param int M: second harmonic term; can be -L <= M <= L |
---|
2343 | :returns True if allowed |
---|
2344 | ''' |
---|
2345 | if not L%2 and abs(M) <= L: |
---|
2346 | if SGLaue == '0': #cylindrical symmetry |
---|
2347 | if M == 0: return True |
---|
2348 | elif SGLaue == '-1': |
---|
2349 | return True |
---|
2350 | elif SGLaue == '2/m': |
---|
2351 | if not abs(M)%2: return True |
---|
2352 | elif SGLaue == 'mmm': |
---|
2353 | if not abs(M)%2 and M >= 0: return True |
---|
2354 | elif SGLaue == '4/m': |
---|
2355 | if not abs(M)%4: return True |
---|
2356 | elif SGLaue == '4/mmm': |
---|
2357 | if not abs(M)%4 and M >= 0: return True |
---|
2358 | elif SGLaue in ['3R','3']: |
---|
2359 | if not abs(M)%3: return True |
---|
2360 | elif SGLaue in ['3mR','3m1','31m']: |
---|
2361 | if not abs(M)%3 and M >= 0: return True |
---|
2362 | elif SGLaue == '6/m': |
---|
2363 | if not abs(M)%6: return True |
---|
2364 | elif SGLaue == '6/mmm': |
---|
2365 | if not abs(M)%6 and M >= 0: return True |
---|
2366 | elif SGLaue in ['m3']: #cubics use different Fourier expansion than those above |
---|
2367 | if M > 0: |
---|
2368 | if L%12 == 2: |
---|
2369 | if M <= L//12: return True |
---|
2370 | else: |
---|
2371 | if M <= L//12+1: return True |
---|
2372 | elif SGLaue in ['m3m']: |
---|
2373 | if M > 0: |
---|
2374 | if L%12 == 2: |
---|
2375 | if M <= L//12: return True |
---|
2376 | else: |
---|
2377 | if M <= L//12+1: return True |
---|
2378 | return False |
---|
2379 | |
---|
2380 | def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True): |
---|
2381 | 'needs doc string' |
---|
2382 | coeffNames = [] |
---|
2383 | for iord in [2*i+2 for i in range(L//2)]: |
---|
2384 | for m in [i-iord for i in range(2*iord+1)]: |
---|
2385 | if OdfChk(SamSym,iord,m): |
---|
2386 | for n in [i-iord for i in range(2*iord+1)]: |
---|
2387 | if OdfChk(SGLaue,iord,n): |
---|
2388 | if IfLMN: |
---|
2389 | coeffNames.append('C(%d,%d,%d)'%(iord,m,n)) |
---|
2390 | else: |
---|
2391 | coeffNames.append('C(%d,%d)'%(iord,n)) |
---|
2392 | return coeffNames |
---|
2393 | |
---|
2394 | def CrsAng(H,cell,SGData): |
---|
2395 | 'needs doc string' |
---|
2396 | a,b,c,al,be,ga = cell |
---|
2397 | SQ3 = 1.732050807569 |
---|
2398 | H1 = np.array([1,0,0]) |
---|
2399 | H2 = np.array([0,1,0]) |
---|
2400 | H3 = np.array([0,0,1]) |
---|
2401 | H4 = np.array([1,1,1]) |
---|
2402 | G,g = cell2Gmat(cell) |
---|
2403 | Laue = SGData['SGLaue'] |
---|
2404 | Naxis = SGData['SGUniq'] |
---|
2405 | if len(H.shape) == 1: |
---|
2406 | DH = np.inner(H,np.inner(G,H)) |
---|
2407 | else: |
---|
2408 | DH = np.array([np.inner(h,np.inner(G,h)) for h in H]) |
---|
2409 | if Laue == '2/m': |
---|
2410 | if Naxis == 'a': |
---|
2411 | DR = np.inner(H1,np.inner(G,H1)) |
---|
2412 | DHR = np.inner(H,np.inner(G,H1)) |
---|
2413 | elif Naxis == 'b': |
---|
2414 | DR = np.inner(H2,np.inner(G,H2)) |
---|
2415 | DHR = np.inner(H,np.inner(G,H2)) |
---|
2416 | else: |
---|
2417 | DR = np.inner(H3,np.inner(G,H3)) |
---|
2418 | DHR = np.inner(H,np.inner(G,H3)) |
---|
2419 | elif Laue in ['R3','R3m']: |
---|
2420 | DR = np.inner(H4,np.inner(G,H4)) |
---|
2421 | DHR = np.inner(H,np.inner(G,H4)) |
---|
2422 | else: |
---|
2423 | DR = np.inner(H3,np.inner(G,H3)) |
---|
2424 | DHR = np.inner(H,np.inner(G,H3)) |
---|
2425 | DHR /= np.sqrt(DR*DH) |
---|
2426 | phi = np.where(DHR <= 1.0,acosd(DHR),0.0) |
---|
2427 | if Laue == '-1': |
---|
2428 | BA = H.T[1]*a/(b-H.T[0]*cosd(ga)) |
---|
2429 | BB = H.T[0]*sind(ga)**2 |
---|
2430 | elif Laue == '2/m': |
---|
2431 | if Naxis == 'a': |
---|
2432 | BA = H.T[2]*b/(c-H.T[1]*cosd(al)) |
---|
2433 | BB = H.T[1]*sind(al)**2 |
---|
2434 | elif Naxis == 'b': |
---|
2435 | BA = H.T[0]*c/(a-H.T[2]*cosd(be)) |
---|
2436 | BB = H.T[2]*sind(be)**2 |
---|
2437 | else: |
---|
2438 | BA = H.T[1]*a/(b-H.T[0]*cosd(ga)) |
---|
2439 | BB = H.T[0]*sind(ga)**2 |
---|
2440 | elif Laue in ['mmm','4/m','4/mmm']: |
---|
2441 | BA = H.T[1]*a |
---|
2442 | BB = H.T[0]*b |
---|
2443 | elif Laue in ['3R','3mR']: |
---|
2444 | BA = H.T[0]+H.T[1]-2.0*H.T[2] |
---|
2445 | BB = SQ3*(H.T[0]-H.T[1]) |
---|
2446 | elif Laue in ['m3','m3m']: |
---|
2447 | BA = H.T[1] |
---|
2448 | BB = H.T[0] |
---|
2449 | else: |
---|
2450 | BA = H.T[0]+2.0*H.T[1] |
---|
2451 | BB = SQ3*H.T[0] |
---|
2452 | beta = atan2d(BA,BB) |
---|
2453 | return phi,beta |
---|
2454 | |
---|
2455 | def SamAng(Tth,Gangls,Sangl,IFCoup): |
---|
2456 | """Compute sample orientation angles vs laboratory coord. system |
---|
2457 | |
---|
2458 | :param Tth: Signed theta |
---|
2459 | :param Gangls: Sample goniometer angles phi,chi,omega,azmuth |
---|
2460 | :param Sangl: Sample angle zeros om-0, chi-0, phi-0 |
---|
2461 | :param IFCoup: True if omega & 2-theta coupled in CW scan |
---|
2462 | :returns: |
---|
2463 | psi,gam: Sample odf angles |
---|
2464 | dPSdA,dGMdA: Angle zero derivatives |
---|
2465 | """ |
---|
2466 | |
---|
2467 | if IFCoup: |
---|
2468 | GSomeg = sind(Gangls[2]+Tth) |
---|
2469 | GComeg = cosd(Gangls[2]+Tth) |
---|
2470 | else: |
---|
2471 | GSomeg = sind(Gangls[2]) |
---|
2472 | GComeg = cosd(Gangls[2]) |
---|
2473 | GSTth = sind(Tth) |
---|
2474 | GCTth = cosd(Tth) |
---|
2475 | GSazm = sind(Gangls[3]) |
---|
2476 | GCazm = cosd(Gangls[3]) |
---|
2477 | GSchi = sind(Gangls[1]) |
---|
2478 | GCchi = cosd(Gangls[1]) |
---|
2479 | GSphi = sind(Gangls[0]+Sangl[2]) |
---|
2480 | GCphi = cosd(Gangls[0]+Sangl[2]) |
---|
2481 | SSomeg = sind(Sangl[0]) |
---|
2482 | SComeg = cosd(Sangl[0]) |
---|
2483 | SSchi = sind(Sangl[1]) |
---|
2484 | SCchi = cosd(Sangl[1]) |
---|
2485 | AT = -GSTth*GComeg+GCTth*GCazm*GSomeg |
---|
2486 | BT = GSTth*GSomeg+GCTth*GCazm*GComeg |
---|
2487 | CT = -GCTth*GSazm*GSchi |
---|
2488 | DT = -GCTth*GSazm*GCchi |
---|
2489 | |
---|
2490 | BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi |
---|
2491 | BC2 = DT-BT*GSchi |
---|
2492 | BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi |
---|
2493 | |
---|
2494 | BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg |
---|
2495 | psi = acosd(BC) |
---|
2496 | |
---|
2497 | BD = 1.0-BC**2 |
---|
2498 | C = np.where(BD>1.e-6,rpd/np.sqrt(BD),0.) |
---|
2499 | dPSdA = [-C*(-BC1*SSomeg*SCchi-BC2*SSomeg*SSchi-BC3*SComeg), |
---|
2500 | -C*(-BC1*SComeg*SSchi+BC2*SComeg*SCchi), |
---|
2501 | -C*(-BC1*SSomeg-BC3*SComeg*SCchi)] |
---|
2502 | |
---|
2503 | BA = -BC1*SSchi+BC2*SCchi |
---|
2504 | BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg |
---|
2505 | gam = atan2d(BB,BA) |
---|
2506 | |
---|
2507 | BD = (BA**2+BB**2)/rpd |
---|
2508 | |
---|
2509 | dBAdO = 0 |
---|
2510 | dBAdC = -BC1*SCchi-BC2*SSchi |
---|
2511 | dBAdF = BC3*SSchi |
---|
2512 | |
---|
2513 | dBBdO = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg |
---|
2514 | dBBdC = -BC1*SSomeg*SSchi+BC2*SSomeg*SCchi |
---|
2515 | dBBdF = BC1*SComeg-BC3*SSomeg*SCchi |
---|
2516 | |
---|
2517 | dGMdA = np.where(BD > 1.e-6,[(BA*dBBdO-BB*dBAdO)/BD,(BA*dBBdC-BB*dBAdC)/BD, \ |
---|
2518 | (BA*dBBdF-BB*dBAdF)/BD],[np.zeros_like(BD),np.zeros_like(BD),np.zeros_like(BD)]) |
---|
2519 | |
---|
2520 | return psi,gam,dPSdA,dGMdA |
---|
2521 | |
---|
2522 | BOH = { |
---|
2523 | 'L=2':[[],[],[]], |
---|
2524 | 'L=4':[[0.30469720,0.36418281],[],[]], |
---|
2525 | 'L=6':[[-0.14104740,0.52775103],[],[]], |
---|
2526 | 'L=8':[[0.28646862,0.21545346,0.32826995],[],[]], |
---|
2527 | 'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]], |
---|
2528 | 'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402], |
---|
2529 | [0.09298802,-0.23773812,0.49446631,0.0],[]], |
---|
2530 | 'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]], |
---|
2531 | 'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495], |
---|
2532 | [0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]], |
---|
2533 | 'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037], |
---|
2534 | [-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]], |
---|
2535 | 'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563], |
---|
2536 | [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]], |
---|
2537 | 'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583], |
---|
2538 | [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]], |
---|
2539 | 'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505], |
---|
2540 | [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0], |
---|
2541 | [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]], |
---|
2542 | 'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594], |
---|
2543 | [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]], |
---|
2544 | 'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501], |
---|
2545 | [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0], |
---|
2546 | [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]], |
---|
2547 | 'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594], |
---|
2548 | [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0], |
---|
2549 | [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]], |
---|
2550 | 'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501], |
---|
2551 | [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0], |
---|
2552 | [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]], |
---|
2553 | 'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508], |
---|
2554 | [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0], |
---|
2555 | [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]] |
---|
2556 | } |
---|
2557 | |
---|
2558 | Lnorm = lambda L: 4.*np.pi/(2.0*L+1.) |
---|
2559 | |
---|
2560 | def GetKcl(L,N,SGLaue,phi,beta): |
---|
2561 | 'needs doc string' |
---|
2562 | import pytexture as ptx |
---|
2563 | if SGLaue in ['m3','m3m']: |
---|
2564 | if 'array' in str(type(phi)) and np.any(phi.shape): |
---|
2565 | Kcl = np.zeros_like(phi) |
---|
2566 | else: |
---|
2567 | Kcl = 0. |
---|
2568 | for j in range(0,L+1,4): |
---|
2569 | im = j//4 |
---|
2570 | if 'array' in str(type(phi)) and np.any(phi.shape): |
---|
2571 | pcrs = ptx.pyplmpsi(L,j,len(phi),phi)[0] |
---|
2572 | else: |
---|
2573 | pcrs = ptx.pyplmpsi(L,j,1,phi)[0] |
---|
2574 | Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta) |
---|
2575 | else: |
---|
2576 | if 'array' in str(type(phi)) and np.any(phi.shape): |
---|
2577 | pcrs = ptx.pyplmpsi(L,N,len(phi),phi)[0] |
---|
2578 | else: |
---|
2579 | pcrs = ptx.pyplmpsi(L,N,1,phi)[0] |
---|
2580 | pcrs *= RSQ2PI |
---|
2581 | if N: |
---|
2582 | pcrs *= SQ2 |
---|
2583 | if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']: |
---|
2584 | if SGLaue in ['3mR','3m1','31m']: |
---|
2585 | if N%6 == 3: |
---|
2586 | Kcl = pcrs*sind(N*beta) |
---|
2587 | else: |
---|
2588 | Kcl = pcrs*cosd(N*beta) |
---|
2589 | else: |
---|
2590 | Kcl = pcrs*cosd(N*beta) |
---|
2591 | else: |
---|
2592 | Kcl = pcrs*(cosd(N*beta)+sind(N*beta)) |
---|
2593 | return Kcl |
---|
2594 | |
---|
2595 | def GetKsl(L,M,SamSym,psi,gam): |
---|
2596 | 'needs doc string' |
---|
2597 | import pytexture as ptx |
---|
2598 | if 'array' in str(type(psi)) and np.any(psi.shape): |
---|
2599 | psrs,dpdps = ptx.pyplmpsi(L,M,len(psi),psi) |
---|
2600 | else: |
---|
2601 | psrs,dpdps = ptx.pyplmpsi(L,M,1,psi) |
---|
2602 | psrs *= RSQ2PI |
---|
2603 | dpdps *= RSQ2PI |
---|
2604 | if M: |
---|
2605 | psrs *= SQ2 |
---|
2606 | dpdps *= SQ2 |
---|
2607 | if SamSym in ['mmm',]: |
---|
2608 | dum = cosd(M*gam) |
---|
2609 | Ksl = psrs*dum |
---|
2610 | dKsdp = dpdps*dum |
---|
2611 | dKsdg = -psrs*M*sind(M*gam) |
---|
2612 | else: |
---|
2613 | dum = cosd(M*gam)+sind(M*gam) |
---|
2614 | Ksl = psrs*dum |
---|
2615 | dKsdp = dpdps*dum |
---|
2616 | dKsdg = psrs*M*(-sind(M*gam)+cosd(M*gam)) |
---|
2617 | return Ksl,dKsdp,dKsdg |
---|
2618 | |
---|
2619 | def GetKclKsl(L,N,SGLaue,psi,phi,beta): |
---|
2620 | """ |
---|
2621 | This is used for spherical harmonics description of preferred orientation; |
---|
2622 | cylindrical symmetry only (M=0) and no sample angle derivatives returned |
---|
2623 | """ |
---|
2624 | import pytexture as ptx |
---|
2625 | Ksl,x = ptx.pyplmpsi(L,0,1,psi) |
---|
2626 | Ksl *= RSQ2PI |
---|
2627 | if SGLaue in ['m3','m3m']: |
---|
2628 | Kcl = 0.0 |
---|
2629 | for j in range(0,L+1,4): |
---|
2630 | im = j//4 |
---|
2631 | pcrs,dum = ptx.pyplmpsi(L,j,1,phi) |
---|
2632 | Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta) |
---|
2633 | else: |
---|
2634 | pcrs,dum = ptx.pyplmpsi(L,N,1,phi) |
---|
2635 | pcrs *= RSQ2PI |
---|
2636 | if N: |
---|
2637 | pcrs *= SQ2 |
---|
2638 | if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']: |
---|
2639 | if SGLaue in ['3mR','3m1','31m']: |
---|
2640 | if N%6 == 3: |
---|
2641 | Kcl = pcrs*sind(N*beta) |
---|
2642 | else: |
---|
2643 | Kcl = pcrs*cosd(N*beta) |
---|
2644 | else: |
---|
2645 | Kcl = pcrs*cosd(N*beta) |
---|
2646 | else: |
---|
2647 | Kcl = pcrs*(cosd(N*beta)+sind(N*beta)) |
---|
2648 | return Kcl*Ksl,Lnorm(L) |
---|
2649 | |
---|
2650 | def H2ThPh(H,Bmat,Q): |
---|
2651 | '''Convert HKL to spherical polar & azimuth angles |
---|
2652 | ''' |
---|
2653 | A,V = G2mth.Q2AVdeg(Q) |
---|
2654 | QR,R = G2mth.make2Quat(V,np.array([0.,0.,1.0])) |
---|
2655 | QA = G2mth.AVdeg2Q(A,np.array([0.,0.,1.0])) |
---|
2656 | Q2 = G2mth.prodQQ(QA,QR) |
---|
2657 | Qmat = G2mth.Q2Mat(Q2) |
---|
2658 | CH = np.inner(H,Bmat) |
---|
2659 | CH = np.inner(CH,Qmat) |
---|
2660 | CH /= nl.norm(CH,axis=2)[:,:,nxs] |
---|
2661 | H3 = np.array([0,0,1.]) |
---|
2662 | DHR = np.inner(CH,H3) |
---|
2663 | Ph = np.where(DHR <= 1.0,acosd(DHR),0.0) #polar angle 0<=Ph<=180. |
---|
2664 | TH = CH*np.array([1.,1.,0.])[nxs,nxs,:] #projection of CH onto xy plane |
---|
2665 | TH /= nl.norm(TH,axis=2)[:,:,nxs] |
---|
2666 | Th = atan2d(TH[:,:,0],TH[:,:,1]) #azimuth angle 0<=Th<360< |
---|
2667 | Th = np.where(Th<0.,Th+360.,Th) |
---|
2668 | return Th,Ph |
---|
2669 | |
---|
2670 | def SHarmcal(SytSym,SHFln,psi,gam): |
---|
2671 | '''Perform a surface spherical harmonics computation. |
---|
2672 | Note that the the number of gam values must either be 1 or must match psi |
---|
2673 | |
---|
2674 | :param str SytSym: sit symmetry - only looking for cubics - remove this |
---|
2675 | :param dict SHFln: spherical harmonics coefficients; key has L & M |
---|
2676 | :param float/array psi: Azimuthal coordinate 0 <= Th <= 360 |
---|
2677 | :param float/array gam: Polar coordinate 0<= Ph <= 180 |
---|
2678 | |
---|
2679 | :returns array SHVal: spherical harmonics array for psi,gam values |
---|
2680 | ''' |
---|
2681 | SHVal = np.ones_like(psi) |
---|
2682 | for term in SHFln: |
---|
2683 | if 'C(' in term[:2]: |
---|
2684 | l,m = eval(term.strip('C').strip('c')) |
---|
2685 | if SytSym in ['m3m','m3','43m','432','23'] or 'c' in term: |
---|
2686 | Ksl = CubicSHarm(l,m,psi,gam) |
---|
2687 | else: |
---|
2688 | p = SHFln[term][2] |
---|
2689 | Ksl = SphHarmAng(l,m,p,psi,gam) |
---|
2690 | SHVal += SHFln[term][0]*Ksl |
---|
2691 | return SHVal |
---|
2692 | |
---|
2693 | def SphHarmAng(L,M,P,Th,Ph): |
---|
2694 | ''' Compute spherical harmonics values using scipy.special.sph_harm |
---|
2695 | |
---|
2696 | :param int L: degree of the harmonic (L >= 0) |
---|
2697 | :param int M: order number (\\|M\\| <= L) |
---|
2698 | :param int P: sign flag = -1 or 1 |
---|
2699 | :param float/array Th: Azimuthal coordinate 0 <= Th <= 360 |
---|
2700 | :param float/array Ph: Polar coordinate 0<= Ph <= 180 |
---|
2701 | |
---|
2702 | :returns ylmp value/array: as reals |
---|
2703 | ''' |
---|
2704 | |
---|
2705 | ylmp = spsp.sph_harm(M,L,rpd*Th,rpd*Ph) #wants radians; order then degree |
---|
2706 | |
---|
2707 | if M == 0: |
---|
2708 | return np.real(ylmp)/SQ2 |
---|
2709 | if P>0: |
---|
2710 | return np.real(ylmp) |
---|
2711 | else: |
---|
2712 | return np.imag(ylmp) |
---|
2713 | |
---|
2714 | def CubicSHarm(L,M,Th,Ph): |
---|
2715 | '''Calculation of the cubic harmonics given in Table 3 in M.Kara & K. Kurki-Suonio, |
---|
2716 | Acta Cryt. A37, 201 (1981). For L = 14,20 only for m3m from F.M. Mueller and M.G. Priestley, |
---|
2717 | Phys Rev 148, 638 (1966) |
---|
2718 | |
---|
2719 | :param int L: degree of the harmonic (L >= 0) |
---|
2720 | :param int M: order number [\\|M\\| <= L] |
---|
2721 | :param float/array Th: Azimuthal coordinate 0 <= Th <= 360 |
---|
2722 | :param float/array Ph: Polar coordinate 0<= Ph <= 180 |
---|
2723 | |
---|
2724 | :returns klm value/array: cubic harmonics |
---|
2725 | |
---|
2726 | ''' |
---|
2727 | if L == 0: |
---|
2728 | return SphHarmAng(L,M,1,Th,Ph) |
---|
2729 | elif L == 3: |
---|
2730 | return SphHarmAng(3,2,-1,Th,Ph) |
---|
2731 | elif L == 4: |
---|
2732 | klm = 0.5*np.sqrt(7.0/3.0)*SphHarmAng(4,0,1,Th,Ph) |
---|
2733 | klm = 0.5*np.sqrt(5.0/3.0)*SphHarmAng(4,4,1,Th,Ph) |
---|
2734 | elif L == 6: |
---|
2735 | if M == 1: |
---|
2736 | klm = 0.5*np.sqrt(0.5)*SphHarmAng(6,0,1,Th,Ph) |
---|
2737 | klm -= 0.5*np.sqrt(7.0/2.0)*SphHarmAng(6,4,1,Th,Ph) |
---|
2738 | else: |
---|
2739 | klm = 0.25*np.sqrt(11.0)*SphHarmAng(6,2,1,Th,Ph) |
---|
2740 | klm -= 0.25*np.sqrt(5.0)*SphHarmAng(6,6,1,Th,Ph) |
---|
2741 | elif L == 7: |
---|
2742 | klm = 0.5*np.sqrt(13./6.)*SphHarmAng(7,2,-1,Th,Ph) |
---|
2743 | klm = 0.5*np.sqrt(11./6.)*SphHarmAng(7,6,-1,Th,Ph) |
---|
2744 | elif L == 8: |
---|
2745 | klm = 0.125*np.sqrt(33.)*SphHarmAng(8,0,1,Th,Ph) |
---|
2746 | klm += 0.25*np.sqrt(7./3.)*SphHarmAng(8,4,1,Th,Ph) |
---|
2747 | klm += 0.125*np.sqrt(65./3.)*SphHarmAng(8,8,1,Th,Ph) |
---|
2748 | elif L == 9: |
---|
2749 | if M == 1: |
---|
2750 | klm = 0.25*np.sqrt(3.)*SphHarmAng(9,2,-1,Th,Ph) |
---|
2751 | klm -= 0.25*np.sqrt(13.)*SphHarmAng(9,6,-1,Th,Ph) |
---|
2752 | else: |
---|
2753 | klm = 0.5*np.sqrt(17./6.)*SphHarmAng(9,4,-1,Th,Ph) |
---|
2754 | klm -= 0.5*np.sqrt(7./6.)*SphHarmAng(9,8,-1,Th,Ph) |
---|
2755 | elif L == 10: |
---|
2756 | if M == 1: |
---|
2757 | klm = 0.125*np.sqrt(65./6.)*SphHarmAng(10,0,1,Th,Ph) |
---|
2758 | klm -= 0.25*np.sqrt(11.0/2.0)*SphHarmAng(10,4,1,Th,Ph) |
---|
2759 | klm -= 0.125*np.sqrt(187.0/6.0)*SphHarmAng(10,8,1,Th,Ph) |
---|
2760 | else: |
---|
2761 | klm = 0.125*np.sqrt(247./6.)*SphHarmAng(10,2,1,Th,Ph) |
---|
2762 | klm += (1./16.)*np.sqrt(19./3.)*SphHarmAng(10,6,1,Th,Ph) |
---|
2763 | klm -= (1./16.)*np.sqrt(85.)*SphHarmAng(10,10,1,Th,Ph) |
---|
2764 | #only m3m cubics from here down; from F.M. Mueller and M.G. Priestley, Phys Rev 148, 638 (1966) |
---|
2765 | elif L == 12: |
---|
2766 | if M == 1: |
---|
2767 | klm = 0.69550265*SphHarmAng(12,0,1,Th,Ph) |
---|
2768 | klm += 0.31412555*SphHarmAng(12,4,1,Th,Ph) |
---|
2769 | klm += 0.34844954*SphHarmAng(12,8,1,Th,Ph) |
---|
2770 | klm += 0.54422797*SphHarmAng(12,12,1,Th,Ph) |
---|
2771 | else: |
---|
2772 | klm = 0.55897937*SphHarmAng(12,4,1,Th,Ph) |
---|
2773 | klm -= 0.80626751*SphHarmAng(12,8,1,Th,Ph) |
---|
2774 | klm += 0.19358400*SphHarmAng(12,12,1,Th,Ph) |
---|
2775 | elif L == 14: |
---|
2776 | klm = 0.44009645*SphHarmAng(14,0,1,Th,Ph) |
---|
2777 | klm -= 0.45768183*SphHarmAng(14,4,1,Th,Ph) |
---|
2778 | klm -= 0.49113230*SphHarmAng(14,8,1,Th,Ph) |
---|
2779 | klm -= 0.59634848*SphHarmAng(14,12,1,Th,Ph) |
---|
2780 | elif L == 16: |
---|
2781 | if M == 1: |
---|
2782 | klm = 0.68136168*SphHarmAng(16,0,1,Th,Ph) |
---|
2783 | klm += 0.27586801*SphHarmAng(16,4,1,Th,Ph) |
---|
2784 | klm += 0.29048987*SphHarmAng(16,8,1,Th,Ph) |
---|
2785 | klm += 0.32756975*SphHarmAng(16,12,1,Th,Ph) |
---|
2786 | klm += 0.51764542*SphHarmAng(16,16,1,Th,Ph) |
---|
2787 | else: |
---|
2788 | klm = 0.63704821*SphHarmAng(16,4,1,Th,Ph) |
---|
2789 | klm -= 0.32999033*SphHarmAng(16,8,1,Th,Ph) |
---|
2790 | klm -= 0.64798073*SphHarmAng(16,12,1,Th,Ph) |
---|
2791 | klm += 0.25572816*SphHarmAng(16,16,1,Th,Ph) |
---|
2792 | elif L == 18: |
---|
2793 | if M == 1: |
---|
2794 | klm = 0.45791513*SphHarmAng(18,0,1,Th,Ph) |
---|
2795 | klm -= 0.38645598*SphHarmAng(18,4,1,Th,Ph) |
---|
2796 | klm -= 0.40209462*SphHarmAng(18,8,1,Th,Ph) |
---|
2797 | klm -= 0.43746593*SphHarmAng(18,12,1,Th,Ph) |
---|
2798 | klm -= 0.53657149*SphHarmAng(18,16,1,Th,Ph) |
---|
2799 | else: |
---|
2800 | klm = 0.14872751*SphHarmAng(18,4,1,Th,Ph) |
---|
2801 | klm -= 0.63774601*SphHarmAng(18,8,1,Th,Ph) |
---|
2802 | klm += 0.72334167*SphHarmAng(18,12,1,Th,Ph) |
---|
2803 | klm -= 0.21894515*SphHarmAng(18,16,1,Th,Ph) |
---|
2804 | elif L == 20: |
---|
2805 | if M == 1: |
---|
2806 | klm = 0.67141495*SphHarmAng(20,0,1,Th,Ph) |
---|
2807 | klm += 0.24982619*SphHarmAng(20,4,1,Th,Ph) |
---|
2808 | klm += 0.25782846*SphHarmAng(20,8,1,Th,Ph) |
---|
2809 | klm += 0.27469333*SphHarmAng(20,12,1,Th,Ph) |
---|
2810 | klm += 0.31248919*SphHarmAng(20,16,1,Th,Ph) |
---|
2811 | klm += 0.49719956*SphHarmAng(20,20,1,Th,Ph) |
---|
2812 | else: |
---|
2813 | klm = 0.66299538*SphHarmAng(20,4,1,Th,Ph) |
---|
2814 | klm -= 0.11295259*SphHarmAng(20,8,1,Th,Ph) |
---|
2815 | klm -= 0.42738441*SphHarmAng(20,12,1,Th,Ph) |
---|
2816 | klm -= 0.52810433*SphHarmAng(20,16,1,Th,Ph) |
---|
2817 | klm += 0.29347435*SphHarmAng(20,20,1,Th,Ph) |
---|
2818 | else: #shouldn't happen |
---|
2819 | return 0.0 |
---|
2820 | return klm |
---|
2821 | |
---|
2822 | def Glnh(SHCoef,psi,gam,SamSym): |
---|
2823 | 'needs doc string' |
---|
2824 | import pytexture as ptx |
---|
2825 | |
---|
2826 | Fln = np.zeros(len(SHCoef)) |
---|
2827 | for i,term in enumerate(SHCoef): |
---|
2828 | l,m,n = eval(term.strip('C')) |
---|
2829 | pcrs,dum = ptx.pyplmpsi(l,m,1,psi) |
---|
2830 | pcrs *= RSQPI |
---|
2831 | if m == 0: |
---|
2832 | pcrs /= SQ2 |
---|
2833 | if SamSym in ['mmm',]: |
---|
2834 | Ksl = pcrs*cosd(m*gam) |
---|
2835 | else: |
---|
2836 | Ksl = pcrs*(cosd(m*gam)+sind(m*gam)) |
---|
2837 | Fln[i] = SHCoef[term]*Ksl*Lnorm(l) |
---|
2838 | ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln)))) |
---|
2839 | return ODFln |
---|
2840 | |
---|
2841 | def Flnh(SHCoef,phi,beta,SGData): |
---|
2842 | 'needs doc string' |
---|
2843 | import pytexture as ptx |
---|
2844 | |
---|
2845 | Fln = np.zeros(len(SHCoef)) |
---|
2846 | for i,term in enumerate(SHCoef): |
---|
2847 | l,m,n = eval(term.strip('C')) |
---|
2848 | if SGData['SGLaue'] in ['m3','m3m']: |
---|
2849 | Kcl = 0.0 |
---|
2850 | for j in range(0,l+1,4): |
---|
2851 | im = j//4 |
---|
2852 | pcrs,dum = ptx.pyplmpsi(l,j,1,phi) |
---|
2853 | Kcl += BOH['L='+str(l)][n-1][im]*pcrs*cosd(j*beta) |
---|
2854 | else: #all but cubic |
---|
2855 | pcrs,dum = ptx.pyplmpsi(l,n,1,phi) |
---|
2856 | pcrs *= RSQPI |
---|
2857 | if n == 0: |
---|
2858 | pcrs /= SQ2 |
---|
2859 | if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']: |
---|
2860 | if SGData['SGLaue'] in ['3mR','3m1','31m']: |
---|
2861 | if n%6 == 3: |
---|
2862 | Kcl = pcrs*sind(n*beta) |
---|
2863 | else: |
---|
2864 | Kcl = pcrs*cosd(n*beta) |
---|
2865 | else: |
---|
2866 | Kcl = pcrs*cosd(n*beta) |
---|
2867 | else: |
---|
2868 | Kcl = pcrs*(cosd(n*beta)+sind(n*beta)) |
---|
2869 | Fln[i] = SHCoef[term]*Kcl*Lnorm(l) |
---|
2870 | ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln)))) |
---|
2871 | return ODFln |
---|
2872 | |
---|
2873 | def polfcal(ODFln,SamSym,psi,gam): |
---|
2874 | '''Perform a pole figure computation. |
---|
2875 | Note that the the number of gam values must either be 1 or must |
---|
2876 | match psi. Updated for numpy 1.8.0 |
---|
2877 | ''' |
---|
2878 | import pytexture as ptx |
---|
2879 | PolVal = np.ones_like(psi) |
---|
2880 | for term in ODFln: |
---|
2881 | if abs(ODFln[term][1]) > 1.e-3: |
---|
2882 | l,m,n = eval(term.strip('C')) |
---|
2883 | psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi) |
---|
2884 | if SamSym in ['-1','2/m']: |
---|
2885 | if m: |
---|
2886 | Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam)) |
---|
2887 | else: |
---|
2888 | Ksl = RSQPI*psrs/SQ2 |
---|
2889 | else: |
---|
2890 | if m: |
---|
2891 | Ksl = RSQPI*psrs*cosd(m*gam) |
---|
2892 | else: |
---|
2893 | Ksl = RSQPI*psrs/SQ2 |
---|
2894 | PolVal += ODFln[term][1]*Ksl |
---|
2895 | return PolVal |
---|
2896 | |
---|
2897 | def invpolfcal(ODFln,SGData,phi,beta): |
---|
2898 | 'needs doc string' |
---|
2899 | import pytexture as ptx |
---|
2900 | |
---|
2901 | invPolVal = np.ones_like(beta) |
---|
2902 | for term in ODFln: |
---|
2903 | if abs(ODFln[term][1]) > 1.e-3: |
---|
2904 | l,m,n = eval(term.strip('C')) |
---|
2905 | if SGData['SGLaue'] in ['m3','m3m']: |
---|
2906 | Kcl = 0.0 |
---|
2907 | for j in range(0,l+1,4): |
---|
2908 | im = j//4 |
---|
2909 | pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi) |
---|
2910 | Kcl += BOH['L=%d'%(l)][n-1][im]*pcrs*cosd(j*beta) |
---|
2911 | else: #all but cubic |
---|
2912 | pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi) |
---|
2913 | pcrs *= RSQPI |
---|
2914 | if n == 0: |
---|
2915 | pcrs /= SQ2 |
---|
2916 | if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']: |
---|
2917 | if SGData['SGLaue'] in ['3mR','3m1','31m']: |
---|
2918 | if n%6 == 3: |
---|
2919 | Kcl = pcrs*sind(n*beta) |
---|
2920 | else: |
---|
2921 | Kcl = pcrs*cosd(n*beta) |
---|
2922 | else: |
---|
2923 | Kcl = pcrs*cosd(n*beta) |
---|
2924 | else: |
---|
2925 | Kcl = pcrs*(cosd(n*beta)+sind(n*beta)) |
---|
2926 | invPolVal += ODFln[term][1]*Kcl |
---|
2927 | return invPolVal |
---|
2928 | |
---|
2929 | |
---|
2930 | def textureIndex(SHCoef): |
---|
2931 | 'needs doc string' |
---|
2932 | Tindx = 1.0 |
---|
2933 | for term in SHCoef: |
---|
2934 | l = eval(term.strip('C'))[0] |
---|
2935 | Tindx += SHCoef[term]**2/(2.0*l+1.) |
---|
2936 | return Tindx |
---|
2937 | |
---|
2938 | UniqueCellByLaue = [ |
---|
2939 | [['m3','m3m'],(0,)], |
---|
2940 | [['3R','3mR'],(0,3)], |
---|
2941 | [['3','3m1','31m','6/m','6/mmm','4/m','4/mmm'],(0,2)], |
---|
2942 | [['mmm'],(0,1,2)], |
---|
2943 | [['2/m'+'a'],(0,1,2,3)], |
---|
2944 | [['2/m'+'b'],(0,1,2,4)], |
---|
2945 | [['2/m'+'c'],(0,1,2,5)], |
---|
2946 | [['-1'],(0,1,2,3,4,5)], |
---|
2947 | ] |
---|
2948 | '''List the unique cell terms by index for each Laue class''' |
---|
2949 | |
---|
2950 | cellAlbl = ('a','b','c', 'alpha', 'beta', 'gamma') |
---|
2951 | 'ASCII labels for a, b, c, alpha, beta, gamma' |
---|
2952 | |
---|
2953 | cellUlbl = ('a','b','c',u'\u03B1',u'\u03B2',u'\u03B3') |
---|
2954 | 'unicode labels for a, b, c, alpha, beta, gamma' |
---|
2955 | |
---|
2956 | # self-test materials follow. |
---|
2957 | selftestlist = [] |
---|
2958 | '''Defines a list of self-tests''' |
---|
2959 | selftestquiet = True |
---|
2960 | def _ReportTest(): |
---|
2961 | 'Report name and doc string of current routine when ``selftestquiet`` is False' |
---|
2962 | if not selftestquiet: |
---|
2963 | import inspect |
---|
2964 | caller = inspect.stack()[1][3] |
---|
2965 | doc = eval(caller).__doc__ |
---|
2966 | if doc is not None: |
---|
2967 | print('testing '+__file__+' with '+caller+' ('+doc+')') |
---|
2968 | else: |
---|
2969 | print('testing '+__file__()+" with "+caller) |
---|
2970 | NeedTestData = True |
---|
2971 | def TestData(): |
---|
2972 | array = np.array |
---|
2973 | global NeedTestData |
---|
2974 | NeedTestData = False |
---|
2975 | global CellTestData |
---|
2976 | # output from uctbx computed on platform darwin on 2010-05-28 |
---|
2977 | CellTestData = [ |
---|
2978 | # cell, g, G, cell*, V, V* |
---|
2979 | [(4, 4, 4, 90, 90, 90), |
---|
2980 | array([[ 1.60000000e+01, 9.79717439e-16, 9.79717439e-16], |
---|
2981 | [ 9.79717439e-16, 1.60000000e+01, 9.79717439e-16], |
---|
2982 | [ 9.79717439e-16, 9.79717439e-16, 1.60000000e+01]]), array([[ 6.25000000e-02, 3.82702125e-18, 3.82702125e-18], |
---|
2983 | [ 3.82702125e-18, 6.25000000e-02, 3.82702125e-18], |
---|
2984 | [ 3.82702125e-18, 3.82702125e-18, 6.25000000e-02]]), (0.25, 0.25, 0.25, 90.0, 90.0, 90.0), 64.0, 0.015625], |
---|
2985 | # cell, g, G, cell*, V, V* |
---|
2986 | [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130), |
---|
2987 | array([[ 16.81 , -13.70423184, 4.48533243], |
---|
2988 | [-13.70423184, 27.04 , -5.6887143 ], |
---|
2989 | [ 4.48533243, -5.6887143 , 39.69 ]]), array([[ 0.10206349, 0.05083339, -0.00424823], |
---|
2990 | [ 0.05083339, 0.06344997, 0.00334956], |
---|
2991 | [-0.00424823, 0.00334956, 0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445], |
---|
2992 | # cell, g, G, cell*, V, V* |
---|
2993 | [(3.5, 3.5, 6, 90, 90, 120), |
---|
2994 | array([[ 1.22500000e+01, -6.12500000e+00, 1.28587914e-15], |
---|
2995 | [ -6.12500000e+00, 1.22500000e+01, 1.28587914e-15], |
---|
2996 | [ 1.28587914e-15, 1.28587914e-15, 3.60000000e+01]]), array([[ 1.08843537e-01, 5.44217687e-02, 3.36690552e-18], |
---|
2997 | [ 5.44217687e-02, 1.08843537e-01, 3.36690552e-18], |
---|
2998 | [ 3.36690552e-18, 3.36690552e-18, 2.77777778e-02]]), (0.32991443953692895, 0.32991443953692895, 0.16666666666666669, 90.0, 90.0, 60.000000000000021), 63.652867178156257, 0.015710211406520427], |
---|
2999 | ] |
---|
3000 | global CoordTestData |
---|
3001 | CoordTestData = [ |
---|
3002 | # cell, ((frac, ortho),...) |
---|
3003 | ((4,4,4,90,90,90,), [ |
---|
3004 | ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)), |
---|
3005 | ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)), |
---|
3006 | ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)), |
---|
3007 | ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)), |
---|
3008 | ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)), |
---|
3009 | ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)), |
---|
3010 | ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)), |
---|
3011 | ]), |
---|
3012 | # cell, ((frac, ortho),...) |
---|
3013 | ((4.1,5.2,6.3,100,80,130,), [ |
---|
3014 | ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)), |
---|
3015 | ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)), |
---|
3016 | ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)), |
---|
3017 | ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)), |
---|
3018 | ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)), |
---|
3019 | ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)), |
---|
3020 | ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)), |
---|
3021 | ]), |
---|
3022 | # cell, ((frac, ortho),...) |
---|
3023 | ((3.5,3.5,6,90,90,120,), [ |
---|
3024 | ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)), |
---|
3025 | ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)), |
---|
3026 | ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)), |
---|
3027 | ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)), |
---|
3028 | ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)), |
---|
3029 | ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)), |
---|
3030 | ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)), |
---|
3031 | ]), |
---|
3032 | ] |
---|
3033 | global LaueTestData #generated by GSAS |
---|
3034 | LaueTestData = { |
---|
3035 | 'R 3 m':[(4.,4.,6.,90.,90.,120.),((1,0,1,6),(1,0,-2,6),(0,0,3,2),(1,1,0,6),(2,0,-1,6),(2,0,2,6), |
---|
3036 | (1,1,3,12),(1,0,4,6),(2,1,1,12),(2,1,-2,12),(3,0,0,6),(1,0,-5,6),(2,0,-4,6),(3,0,-3,6),(3,0,3,6), |
---|
3037 | (0,0,6,2),(2,2,0,6),(2,1,4,12),(2,0,5,6),(3,1,-1,12),(3,1,2,12),(1,1,6,12),(2,2,3,12),(2,1,-5,12))], |
---|
3038 | 'R 3':[(4.,4.,6.,90.,90.,120.),((1,0,1,6),(1,0,-2,6),(0,0,3,2),(1,1,0,6),(2,0,-1,6),(2,0,2,6),(1,1,3,6), |
---|
3039 | (1,1,-3,6),(1,0,4,6),(3,-1,1,6),(2,1,1,6),(3,-1,-2,6),(2,1,-2,6),(3,0,0,6),(1,0,-5,6),(2,0,-4,6), |
---|
3040 | (2,2,0,6),(3,0,3,6),(3,0,-3,6),(0,0,6,2),(3,-1,4,6),(2,0,5,6),(2,1,4,6),(4,-1,-1,6),(3,1,-1,6), |
---|
3041 | (3,1,2,6),(4,-1,2,6),(2,2,-3,6),(1,1,-6,6),(1,1,6,6),(2,2,3,6),(2,1,-5,6),(3,-1,-5,6))], |
---|
3042 | 'P 3':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(1,0,1,6),(0,0,2,2),(1,0,-1,6),(1,0,2,6),(1,0,-2,6), |
---|
3043 | (1,1,0,6),(0,0,3,2),(1,1,1,6),(1,1,-1,6),(1,0,3,6),(1,0,-3,6),(2,0,0,6),(2,0,-1,6),(1,1,-2,6), |
---|
3044 | (1,1,2,6),(2,0,1,6),(2,0,-2,6),(2,0,2,6),(0,0,4,2),(1,1,-3,6),(1,1,3,6),(1,0,-4,6),(1,0,4,6), |
---|
3045 | (2,0,-3,6),(2,1,0,6),(2,0,3,6),(3,-1,0,6),(2,1,1,6),(3,-1,-1,6),(2,1,-1,6),(3,-1,1,6),(1,1,4,6), |
---|
3046 | (3,-1,2,6),(3,-1,-2,6),(1,1,-4,6),(0,0,5,2),(2,1,2,6),(2,1,-2,6),(3,0,0,6),(3,0,1,6),(2,0,4,6), |
---|
3047 | (2,0,-4,6),(3,0,-1,6),(1,0,-5,6),(1,0,5,6),(3,-1,-3,6),(2,1,-3,6),(2,1,3,6),(3,-1,3,6),(3,0,-2,6), |
---|
3048 | (3,0,2,6),(1,1,5,6),(1,1,-5,6),(2,2,0,6),(3,0,3,6),(3,0,-3,6),(0,0,6,2),(2,0,-5,6),(2,1,-4,6), |
---|
3049 | (2,2,-1,6),(3,-1,-4,6),(2,2,1,6),(3,-1,4,6),(2,1,4,6),(2,0,5,6),(1,0,-6,6),(1,0,6,6),(4,-1,0,6), |
---|
3050 | (3,1,0,6),(3,1,-1,6),(3,1,1,6),(4,-1,-1,6),(2,2,2,6),(4,-1,1,6),(2,2,-2,6),(3,1,2,6),(3,1,-2,6), |
---|
3051 | (3,0,4,6),(3,0,-4,6),(4,-1,-2,6),(4,-1,2,6),(2,2,-3,6),(1,1,6,6),(1,1,-6,6),(2,2,3,6),(3,-1,5,6), |
---|
3052 | (2,1,5,6),(2,1,-5,6),(3,-1,-5,6))], |
---|
3053 | 'P 3 m 1':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(1,0,-1,6),(1,0,1,6),(0,0,2,2),(1,0,-2,6), |
---|
3054 | (1,0,2,6),(1,1,0,6),(0,0,3,2),(1,1,1,12),(1,0,-3,6),(1,0,3,6),(2,0,0,6),(1,1,2,12),(2,0,1,6), |
---|
3055 | (2,0,-1,6),(0,0,4,2),(2,0,-2,6),(2,0,2,6),(1,1,3,12),(1,0,-4,6),(1,0,4,6),(2,0,3,6),(2,1,0,12), |
---|
3056 | (2,0,-3,6),(2,1,1,12),(2,1,-1,12),(1,1,4,12),(2,1,2,12),(0,0,5,2),(2,1,-2,12),(3,0,0,6),(1,0,-5,6), |
---|
3057 | (3,0,1,6),(3,0,-1,6),(1,0,5,6),(2,0,4,6),(2,0,-4,6),(2,1,3,12),(2,1,-3,12),(3,0,-2,6),(3,0,2,6), |
---|
3058 | (1,1,5,12),(3,0,-3,6),(0,0,6,2),(2,2,0,6),(3,0,3,6),(2,1,4,12),(2,2,1,12),(2,0,5,6),(2,1,-4,12), |
---|
3059 | (2,0,-5,6),(1,0,-6,6),(1,0,6,6),(3,1,0,12),(3,1,-1,12),(3,1,1,12),(2,2,2,12),(3,1,2,12), |
---|
3060 | (3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))], |
---|
3061 | 'P 3 1 m':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(0,0,2,2),(1,0,1,12),(1,0,2,12),(1,1,0,6), |
---|
3062 | (0,0,3,2),(1,1,-1,6),(1,1,1,6),(1,0,3,12),(2,0,0,6),(2,0,1,12),(1,1,2,6),(1,1,-2,6),(2,0,2,12), |
---|
3063 | (0,0,4,2),(1,1,-3,6),(1,1,3,6),(1,0,4,12),(2,1,0,12),(2,0,3,12),(2,1,1,12),(2,1,-1,12),(1,1,-4,6), |
---|
3064 | (1,1,4,6),(0,0,5,2),(2,1,-2,12),(2,1,2,12),(3,0,0,6),(1,0,5,12),(2,0,4,12),(3,0,1,12),(2,1,-3,12), |
---|
3065 | (2,1,3,12),(3,0,2,12),(1,1,5,6),(1,1,-5,6),(3,0,3,12),(0,0,6,2),(2,2,0,6),(2,1,-4,12),(2,0,5,12), |
---|
3066 | (2,2,-1,6),(2,2,1,6),(2,1,4,12),(3,1,0,12),(1,0,6,12),(2,2,2,6),(3,1,-1,12),(2,2,-2,6),(3,1,1,12), |
---|
3067 | (3,1,-2,12),(3,0,4,12),(3,1,2,12),(1,1,-6,6),(2,2,3,6),(2,2,-3,6),(1,1,6,6))], |
---|
3068 | } |
---|
3069 | |
---|
3070 | global FLnhTestData |
---|
3071 | FLnhTestData = [{ |
---|
3072 | 'C(4,0,0)': (0.965, 0.42760447), |
---|
3073 | 'C(2,0,0)': (1.0122, -0.80233610), |
---|
3074 | 'C(2,0,2)': (0.0061, 8.37491546E-03), |
---|
3075 | 'C(6,0,4)': (-0.0898, 4.37985696E-02), |
---|
3076 | 'C(6,0,6)': (-0.1369, -9.04081762E-02), |
---|
3077 | 'C(6,0,0)': (0.5935, -0.18234928), |
---|
3078 | 'C(4,0,4)': (0.1872, 0.16358127), |
---|
3079 | 'C(6,0,2)': (0.6193, 0.27573633), |
---|
3080 | 'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]] |
---|
3081 | def test0(): |
---|
3082 | if NeedTestData: TestData() |
---|
3083 | msg = 'test cell2Gmat, fillgmat, Gmat2cell' |
---|
3084 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
3085 | G, g = cell2Gmat(cell) |
---|
3086 | assert np.allclose(G,tG),msg |
---|
3087 | assert np.allclose(g,tg),msg |
---|
3088 | tcell = Gmat2cell(g) |
---|
3089 | assert np.allclose(cell,tcell),msg |
---|
3090 | tcell = Gmat2cell(G) |
---|
3091 | assert np.allclose(tcell,trcell),msg |
---|
3092 | if __name__ == '__main__': selftestlist.append(test0) |
---|
3093 | |
---|
3094 | def test1(): |
---|
3095 | 'test cell2A and A2Gmat' |
---|
3096 | _ReportTest() |
---|
3097 | if NeedTestData: TestData() |
---|
3098 | msg = 'test cell2A and A2Gmat' |
---|
3099 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
3100 | G, g = A2Gmat(cell2A(cell)) |
---|
3101 | assert np.allclose(G,tG),msg |
---|
3102 | assert np.allclose(g,tg),msg |
---|
3103 | if __name__ == '__main__': selftestlist.append(test1) |
---|
3104 | |
---|
3105 | def test2(): |
---|
3106 | 'test Gmat2A, A2cell, A2Gmat, Gmat2cell' |
---|
3107 | _ReportTest() |
---|
3108 | if NeedTestData: TestData() |
---|
3109 | msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell' |
---|
3110 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
3111 | G, g = cell2Gmat(cell) |
---|
3112 | tcell = A2cell(Gmat2A(G)) |
---|
3113 | assert np.allclose(cell,tcell),msg |
---|
3114 | if __name__ == '__main__': selftestlist.append(test2) |
---|
3115 | |
---|
3116 | def test3(): |
---|
3117 | 'test invcell2Gmat' |
---|
3118 | _ReportTest() |
---|
3119 | if NeedTestData: TestData() |
---|
3120 | msg = 'test invcell2Gmat' |
---|
3121 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
3122 | G, g = invcell2Gmat(trcell) |
---|
3123 | assert np.allclose(G,tG),msg |
---|
3124 | assert np.allclose(g,tg),msg |
---|
3125 | if __name__ == '__main__': selftestlist.append(test3) |
---|
3126 | |
---|
3127 | def test4(): |
---|
3128 | 'test calc_rVsq, calc_rV, calc_V' |
---|
3129 | _ReportTest() |
---|
3130 | if NeedTestData: TestData() |
---|
3131 | msg = 'test calc_rVsq, calc_rV, calc_V' |
---|
3132 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
3133 | assert np.allclose(calc_rV(cell2A(cell)),trV), msg |
---|
3134 | assert np.allclose(calc_V(cell2A(cell)),tV), msg |
---|
3135 | if __name__ == '__main__': selftestlist.append(test4) |
---|
3136 | |
---|
3137 | def test5(): |
---|
3138 | 'test A2invcell' |
---|
3139 | _ReportTest() |
---|
3140 | if NeedTestData: TestData() |
---|
3141 | msg = 'test A2invcell' |
---|
3142 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
3143 | rcell = A2invcell(cell2A(cell)) |
---|
3144 | assert np.allclose(rcell,trcell),msg |
---|
3145 | if __name__ == '__main__': selftestlist.append(test5) |
---|
3146 | |
---|
3147 | def test6(): |
---|
3148 | 'test cell2AB' |
---|
3149 | _ReportTest() |
---|
3150 | if NeedTestData: TestData() |
---|
3151 | msg = 'test cell2AB' |
---|
3152 | for (cell,coordlist) in CoordTestData: |
---|
3153 | A,B = cell2AB(cell) |
---|
3154 | for (frac,ortho) in coordlist: |
---|
3155 | to = np.inner(A,frac) |
---|
3156 | tf = np.inner(B,to) |
---|
3157 | assert np.allclose(ortho,to), msg |
---|
3158 | assert np.allclose(frac,tf), msg |
---|
3159 | to = np.sum(A*frac,axis=1) |
---|
3160 | tf = np.sum(B*to,axis=1) |
---|
3161 | assert np.allclose(ortho,to), msg |
---|
3162 | assert np.allclose(frac,tf), msg |
---|
3163 | if __name__ == '__main__': selftestlist.append(test6) |
---|
3164 | |
---|
3165 | def test7(): |
---|
3166 | 'test GetBraviasNum(...) and GenHBravais(...)' |
---|
3167 | _ReportTest() |
---|
3168 | import os.path |
---|
3169 | import sys |
---|
3170 | import GSASIIspc as spc |
---|
3171 | testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp') |
---|
3172 | if os.path.exists(testdir): |
---|
3173 | if testdir not in sys.path: sys.path.insert(0,testdir) |
---|
3174 | import sgtbxlattinp |
---|
3175 | derror = 1e-4 |
---|
3176 | def indexmatch(hklin, hkllist, system): |
---|
3177 | for hklref in hkllist: |
---|
3178 | hklref = list(hklref) |
---|
3179 | # these permutations are far from complete, but are sufficient to |
---|
3180 | # allow the test to complete |
---|
3181 | if system == 'cubic': |
---|
3182 | permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),] |
---|
3183 | elif system == 'monoclinic': |
---|
3184 | permlist = [(1,2,3),(-1,2,-3)] |
---|
3185 | else: |
---|
3186 | permlist = [(1,2,3)] |
---|
3187 | |
---|
3188 | for perm in permlist: |
---|
3189 | hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm] |
---|
3190 | if hkl == hklref: return True |
---|
3191 | if [-i for i in hkl] == hklref: return True |
---|
3192 | else: |
---|
3193 | return False |
---|
3194 | |
---|
3195 | for key in sgtbxlattinp.sgtbx7: |
---|
3196 | spdict = spc.SpcGroup(key) |
---|
3197 | cell = sgtbxlattinp.sgtbx7[key][0] |
---|
3198 | system = spdict[1]['SGSys'] |
---|
3199 | center = spdict[1]['SGLatt'] |
---|
3200 | |
---|
3201 | bravcode = GetBraviasNum(center, system) |
---|
3202 | |
---|
3203 | g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell)) |
---|
3204 | |
---|
3205 | assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key |
---|
3206 | for h,k,l,d,num in g2list: |
---|
3207 | for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]: |
---|
3208 | if abs(d-dref) < derror: |
---|
3209 | if indexmatch((h,k,l,), hkllist, system): |
---|
3210 | break |
---|
3211 | else: |
---|
3212 | assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key) |
---|
3213 | if __name__ == '__main__': selftestlist.append(test7) |
---|
3214 | |
---|
3215 | def test8(): |
---|
3216 | 'test GenHLaue' |
---|
3217 | _ReportTest() |
---|
3218 | import GSASIIspc as spc |
---|
3219 | import sgtbxlattinp |
---|
3220 | derror = 1e-4 |
---|
3221 | dmin = sgtbxlattinp.dmin |
---|
3222 | |
---|
3223 | def indexmatch(hklin, hklref, system, axis): |
---|
3224 | # these permutations are far from complete, but are sufficient to |
---|
3225 | # allow the test to complete |
---|
3226 | if system == 'cubic': |
---|
3227 | permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),] |
---|
3228 | elif system == 'monoclinic' and axis=='b': |
---|
3229 | permlist = [(1,2,3),(-1,2,-3)] |
---|
3230 | elif system == 'monoclinic' and axis=='a': |
---|
3231 | permlist = [(1,2,3),(1,-2,-3)] |
---|
3232 | elif system == 'monoclinic' and axis=='c': |
---|
3233 | permlist = [(1,2,3),(-1,-2,3)] |
---|
3234 | elif system == 'trigonal': |
---|
3235 | permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)] |
---|
3236 | elif system == 'rhombohedral': |
---|
3237 | permlist = [(1,2,3),(2,3,1),(3,1,2)] |
---|
3238 | else: |
---|
3239 | permlist = [(1,2,3)] |
---|
3240 | |
---|
3241 | hklref = list(hklref) |
---|
3242 | for perm in permlist: |
---|
3243 | hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm] |
---|
3244 | if hkl == hklref: return True |
---|
3245 | if [-i for i in hkl] == hklref: return True |
---|
3246 | return False |
---|
3247 | |
---|
3248 | for key in sgtbxlattinp.sgtbx8: |
---|
3249 | spdict = spc.SpcGroup(key)[1] |
---|
3250 | cell = sgtbxlattinp.sgtbx8[key][0] |
---|
3251 | Axis = spdict['SGUniq'] |
---|
3252 | system = spdict['SGSys'] |
---|
3253 | |
---|
3254 | g2list = GenHLaue(dmin,spdict,cell2A(cell)) |
---|
3255 | #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]): |
---|
3256 | # print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1]) |
---|
3257 | # print 'GSAS-II:' |
---|
3258 | # for h,k,l,d in g2list: print ' ',(h,k,l),d |
---|
3259 | # print 'SGTBX:' |
---|
3260 | # for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print ' ',hkllist,dref |
---|
3261 | assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), ( |
---|
3262 | 'Reflection lists differ for %s' % key |
---|
3263 | ) |
---|
3264 | #match = True |
---|
3265 | for h,k,l,d in g2list: |
---|
3266 | for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: |
---|
3267 | if abs(d-dref) < derror: |
---|
3268 | if indexmatch((h,k,l,), hkllist, system, Axis): break |
---|
3269 | else: |
---|
3270 | assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key) |
---|
3271 | #match = False |
---|
3272 | #if not match: |
---|
3273 | #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print ' ',hkllist,dref |
---|
3274 | #print center, Laue, Axis, system |
---|
3275 | if __name__ == '__main__': selftestlist.append(test8) |
---|
3276 | |
---|
3277 | def test9(): |
---|
3278 | 'test GenHLaue' |
---|
3279 | _ReportTest() |
---|
3280 | import GSASIIspc as G2spc |
---|
3281 | if NeedTestData: TestData() |
---|
3282 | for spc in LaueTestData: |
---|
3283 | data = LaueTestData[spc] |
---|
3284 | cell = data[0] |
---|
3285 | hklm = np.array(data[1]) |
---|
3286 | H = hklm[-1][:3] |
---|
3287 | hklO = hklm.T[:3].T |
---|
3288 | A = cell2A(cell) |
---|
3289 | dmin = 1./np.sqrt(calc_rDsq(H,A)) |
---|
3290 | SGData = G2spc.SpcGroup(spc)[1] |
---|
3291 | hkls = np.array(GenHLaue(dmin,SGData,A)) |
---|
3292 | hklN = hkls.T[:3].T |
---|
3293 | #print spc,hklO.shape,hklN.shape |
---|
3294 | err = True |
---|
3295 | for H in hklO: |
---|
3296 | if H not in hklN: |
---|
3297 | print ('%d %s'%(H,' missing from hkl from GSASII')) |
---|
3298 | err = False |
---|
3299 | assert(err) |
---|
3300 | if __name__ == '__main__': selftestlist.append(test9) |
---|
3301 | |
---|
3302 | |
---|
3303 | |
---|
3304 | |
---|
3305 | if __name__ == '__main__': |
---|
3306 | # run self-tests |
---|
3307 | selftestquiet = False |
---|
3308 | for test in selftestlist: |
---|
3309 | test() |
---|
3310 | print ("OK") |
---|