source: trunk/GSASIIspc.py @ 3619

Last change on this file since 3619 was 3619, checked in by vondreele, 4 years ago

clear magData when a selection is made
add row select --> mag space group display in magCells table
fix R -3 m spin selection problem

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 169.9 KB
Line 
1# -*- coding: utf-8 -*-
2"""
3*GSASIIspc: Space group module*
4-------------------------------
5
6Space group interpretation routines. Note that space group information is
7stored in a :ref:`Space Group (SGData)<SGData_table>` object.
8
9"""
10########### SVN repository information ###################
11# $Date: 2018-09-25 16:07:53 +0000 (Tue, 25 Sep 2018) $
12# $Author: vondreele $
13# $Revision: 3619 $
14# $URL: trunk/GSASIIspc.py $
15# $Id: GSASIIspc.py 3619 2018-09-25 16:07:53Z vondreele $
16########### SVN repository information ###################
17from __future__ import division, print_function
18import numpy as np
19import numpy.linalg as nl
20import scipy.optimize as so
21import sys
22import copy
23import os.path as ospath
24
25import GSASIIpath
26GSASIIpath.SetVersionNumber("$Revision: 3619 $")
27
28npsind = lambda x: np.sin(x*np.pi/180.)
29npcosd = lambda x: np.cos(x*np.pi/180.)
30nxs = np.newaxis
31twopi = 2.0*np.pi
32DEBUG = False
33   
34################################################################################
35#### Space group codes
36################################################################################
37
38def SpcGroup(SGSymbol):
39    """
40    Determines cell and symmetry information from a short H-M space group name
41
42    :param SGSymbol: space group symbol (string) with spaces between axial fields
43    :returns: (SGError,SGData)
44   
45       * SGError = 0 for no errors; >0 for errors (see SGErrors below for details)
46       * SGData - is a dict (see :ref:`Space Group object<SGData_table>`) with entries:
47       
48             * 'SpGrp': space group symbol, slightly cleaned up
49             * 'SGFixed': True if space group data can not be changed, e.g. from magnetic cif; otherwise False
50             * 'SGGray': True if 1' in symbol - gray group for mag. incommensurate phases
51             * 'SGLaue':  one of '-1', '2/m', 'mmm', '4/m', '4/mmm', '3R',
52               '3mR', '3', '3m1', '31m', '6/m', '6/mmm', 'm3', 'm3m'
53             * 'SGInv': boolean; True if centrosymmetric, False if not
54             * 'SGLatt': one of 'P', 'A', 'B', 'C', 'I', 'F', 'R'
55             * 'SGUniq': one of 'a', 'b', 'c' if monoclinic, '' otherwise
56             * 'SGCen': cell centering vectors [0,0,0] at least
57             * 'SGOps': symmetry operations as [M,T] so that M*x+T = x'
58             * 'SGSys': one of 'triclinic', 'monoclinic', 'orthorhombic',
59               'tetragonal', 'rhombohedral', 'trigonal', 'hexagonal', 'cubic'
60             * 'SGPolax': one of ' ', 'x', 'y', 'x y', 'z', 'x z', 'y z',
61               'xyz', '111' for arbitrary axes
62             * 'SGPtGrp': one of 32 point group symbols (with some permutations), which
63                is filled by SGPtGroup, is external (KE) part of supersymmetry point group
64             * 'SSGKl': default internal (Kl) part of supersymmetry point group; modified
65                in supersymmetry stuff depending on chosen modulation vector for Mono & Ortho
66             * 'BNSlattsym': BNS lattice symbol & cenering op - used for magnetic structures
67
68    """
69    LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
70    LattSym = ('P','A','B','C','I','F','R')
71    UniqSym = ('','','a','b','c','',)
72    SysSym = ('triclinic','monoclinic','orthorhombic','tetragonal','rhombohedral','trigonal','hexagonal','cubic')
73    SGData = {}
74    if len(SGSymbol.split()) < 2:
75        return SGErrors(0),SGData
76    if ':R' in SGSymbol:
77        SGSymbol = SGSymbol.replace(':',' ')    #get rid of ':' in R space group symbols from some cif files
78    SGData['SGGray'] = False
79    if "1'" in SGSymbol:        #set for incommensurate magnetic
80        SGData['SGGray'] = True
81        SGSymbol = SGSymbol.replace("1'",'')
82    SGSymbol = SGSymbol.split(':')[0]   #remove :1/2 setting symbol from some cif files
83    if '-2' in SGSymbol:    #replace bad but legal symbols with correct equivalents
84        SGSymbol = SGSymbol.replace('-2','m')
85    if SGSymbol.split()[1] =='3/m':
86        SGSymbol = SGSymbol.replace('3/m','-6')
87    import pyspg
88    SGInfo = pyspg.sgforpy(SGSymbol)
89    SGData['SpGrp'] = SGSymbol.strip().lower().capitalize()
90    SGData['SGLaue'] = LaueSym[SGInfo[0]-1]
91    SGData['SGInv'] = bool(SGInfo[1])
92    SGData['SGLatt'] = LattSym[SGInfo[2]-1]
93    SGData['SGUniq'] = UniqSym[SGInfo[3]+1]
94    SGData['SGFixed'] = False
95    SGData['SGOps'] = []
96    SGData['SGGen'] = []
97    for i in range(SGInfo[5]):
98        Mat = np.array(SGInfo[6][i])
99        Trns = np.array(SGInfo[7][i])
100        SGData['SGOps'].append([Mat,Trns])
101        if 'array' in str(type(SGInfo[8])):        #patch for old fortran bin?
102            SGData['SGGen'].append(int(SGInfo[8][i]))
103    SGData['BNSlattsym'] = [LattSym[SGInfo[2]-1],[0,0,0]]
104    lattSpin = []
105    if SGData['SGLatt'] == 'P':
106        SGData['SGCen'] = np.array(([0,0,0],))
107    elif SGData['SGLatt'] == 'A':
108        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5]))
109        lattSpin += [1,]
110    elif SGData['SGLatt'] == 'B':
111        SGData['SGCen'] = np.array(([0,0,0],[.5,0,.5]))
112        lattSpin += [1,]
113    elif SGData['SGLatt'] == 'C':
114        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,0,]))
115        lattSpin += [1,]
116    elif SGData['SGLatt'] == 'I':
117        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,.5]))
118        lattSpin += [1,]
119    elif SGData['SGLatt'] == 'F':
120        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5],[.5,0,.5],[.5,.5,0,]))
121        lattSpin += [1,1,1,1]
122    elif SGData['SGLatt'] == 'R':
123        SGData['SGCen'] = np.array(([0,0,0],[2./3,1./3,1./3],[1./3,2./3,2./3]))
124
125    if SGData['SGInv']:
126        if SGData['SGLaue'] in ['-1','2/m','mmm']:
127            Ibar = 7
128        elif SGData['SGLaue'] in ['4/m','4/mmm']:
129            Ibar = 1
130        elif SGData['SGLaue'] in ['3R','3mR','3','3m1','31m','6/m','6/mmm']:
131            Ibar = 15 #8+4+2+1
132        else:
133            Ibar = 4
134        Ibarx = Ibar&14
135    else:
136        Ibarx = 8
137        if SGData['SGLaue'] in ['-1','2/m','mmm','m3','m3m']:
138            Ibarx = 0
139    moregen = []
140    for i,gen in enumerate(SGData['SGGen']):
141        if SGData['SGLaue'] in ['m3','m3m']:
142            if gen in [1,2,4]:
143                SGData['SGGen'][i] = 4
144            elif gen < 7:
145                SGData['SGGen'][i] = 0
146        elif SGData['SGLaue'] in ['4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm']:
147            if gen == 2:
148                SGData['SGGen'][i] = 4
149            elif gen in [3,5]:
150                SGData['SGGen'][i] = 3
151            elif gen == 6:
152                if SGData['SGLaue'] in ['4/m','4/mmm']:
153                    SGData['SGGen'][i] = 128
154                else:
155                    SGData['SGGen'][i] = 16
156            elif not SGData['SGInv'] and gen == 12:
157                SGData['SGGen'][i] = 8
158            elif (not SGData['SGInv']) and (SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']) and (gen == 1):
159                SGData['SGGen'][i] = 24
160        gen = SGData['SGGen'][i]
161        if gen == 99:
162            gen = 8
163            if SGData['SGLaue'] in ['3m1','31m','6/m','6/mmm']:
164                gen = 3
165            elif SGData['SGLaue'] == 'm3m':
166                gen = 12
167            SGData['SGGen'][i] = gen
168        elif gen == 98:
169            gen = 8
170            if SGData['SGLaue'] in ['3m1','31m','6/m','6/mmm']:
171                gen = 4
172            SGData['SGGen'][i] = gen
173        elif not SGData['SGInv'] and gen in [23,] and SGData['SGLaue'] in ['m3','m3m']:
174            SGData['SGGen'][i] = 24
175        elif gen >= 16 and gen != 128:
176            if not SGData['SGInv']:
177                gen = 31
178            else:
179                gen ^= Ibarx
180            SGData['SGGen'][i] = gen
181        if SGData['SGInv']:
182            if gen < 128:
183                moregen.append(SGData['SGGen'][i]^Ibar)
184            else:
185                moregen.append(1)
186    SGData['SGGen'] += moregen
187    if SGData['SGLaue'] in '-1':
188        SGData['SGSys'] = SysSym[0]
189    elif SGData['SGLaue'] in '2/m':
190        SGData['SGSys'] = SysSym[1]
191    elif SGData['SGLaue'] in 'mmm':
192        SGData['SGSys'] = SysSym[2]
193    elif SGData['SGLaue'] in ['4/m','4/mmm']:
194        SGData['SGSys'] = SysSym[3]
195    elif SGData['SGLaue'] in ['3R','3mR']:
196        SGData['SGSys'] = SysSym[4]
197    elif SGData['SGLaue'] in ['3','3m1','31m']:
198        SGData['SGSys'] = SysSym[5]
199    elif SGData['SGLaue'] in ['6/m','6/mmm']:
200        SGData['SGSys'] = SysSym[6]
201    elif SGData['SGLaue'] in ['m3','m3m']:
202        SGData['SGSys'] = SysSym[7]
203    SGData['SGPolax'] = SGpolar(SGData)
204    SGData['SGPtGrp'],SGData['SSGKl'] = SGPtGroup(SGData)
205
206    if SGData['SGLatt'] == 'R':
207        if SGData['SGPtGrp'] in ['3',]:
208            SGData['SGSpin'] = 3*[1,]
209        elif SGData['SGPtGrp'] in  ['-3','32','3m']:
210            SGData['SGSpin'] = 4*[1,]
211        elif SGData['SGPtGrp'] in  ['-3m',]:
212            SGData['SGSpin'] = 5*[1,]
213       
214    else:
215        if SGData['SGPtGrp'] in ['1','3','23',]:
216            SGData['SGSpin'] = lattSpin+[1,]
217        elif SGData['SGPtGrp'] in ['-1','2','m','4','-4','-3','312','321','3m1','31m','6','-6','432','-43m']:
218            SGData['SGSpin'] = lattSpin+[1,1,]
219        elif SGData['SGPtGrp'] in ['2/m','4/m','422','4mm','-42m','-4m2','-3m1','-31m',
220            '6/m','622','6mm','-6m2','-62m','m3','m3m']:
221            SGData['SGSpin'] = lattSpin+[1,1,1,]
222        else: #'222'-'mmm','4/mmm','6/mmm'
223            SGData['SGSpin'] = lattSpin+[1,1,1,1,]
224    return SGInfo[-1],SGData
225
226def SGErrors(IErr):
227    '''
228    Interprets the error message code from SpcGroup. Used in SpaceGroup.
229   
230    :param IErr: see SGError in :func:`SpcGroup`
231    :returns:
232        ErrString - a string with the error message or "Unknown error"
233    '''
234
235    ErrString = [' ',
236        'Less than 2 operator fields were found',
237        'Illegal Lattice type, not P, A, B, C, I, F or R',
238        'Rhombohedral lattice requires a 3-axis',
239        'Minus sign does not preceed 1, 2, 3, 4 or 6',
240        'Either a 5-axis anywhere or a 3-axis in field not allowed',
241        ' ',
242        'I for COMPUTED GO TO out of range.',
243        'An a-glide mirror normal to A not allowed',
244        'A b-glide mirror normal to B not allowed',
245        'A c-glide mirror normal to C not allowed',
246        'D-glide in a primitive lattice not allowed',
247        'A 4-axis not allowed in the 2nd operator field',
248        'A 6-axis not allowed in the 2nd operator field',
249        'More than 24 matrices needed to define group',
250        ' ',
251        'Improper construction of a rotation operator',
252        'Mirror following a / not allowed',
253        'A translation conflict between operators',
254        'The 2bar operator is not allowed',
255        '3 fields are legal only in R & m3 cubic groups',
256        'Syntax error. Expected I -4 3 d at this point',
257        ' ',
258        'A or B centered tetragonal not allowed',
259        ' ','unknown error in sgroup',' ',' ',' ',
260        'Illegal character in the space group symbol',
261        ]
262    try:
263        return ErrString[IErr]
264    except:
265        return "Unknown error"
266
267def SGpolar(SGData):
268    '''
269    Determine identity of polar axes if any
270    '''
271    POL = ('','x','y','x y','z','x z','y z','xyz','111')
272    NP = [1,2,4]
273    NPZ = [0,1]
274    for M,T in SGData['SGOps']:
275        for i in range(3):
276            if M[i][i] <= 0.: NP[i] = 0
277        if M[0][2] > 0: NPZ[0] = 8
278        if M[1][2] > 0: NPZ[1] = 0
279    NPol = (NP[0]+NP[1]+NP[2]+NPZ[0]*NPZ[1])*(1-int(SGData['SGInv']))
280    return POL[NPol]
281   
282def SGPtGroup(SGData):
283    '''
284    Determine point group of the space group - done after space group symbol has
285    been evaluated by SpcGroup. Only short symbols are allowed
286   
287    :param SGData: from :func SpcGroup
288    :returns: SSGPtGrp & SSGKl (only defaults for Mono & Ortho)
289    '''
290    Flds = SGData['SpGrp'].split()
291    if len(Flds) < 2:
292        return '',[]
293    if SGData['SGLaue'] == '-1':    #triclinic
294        if '-' in Flds[1]:
295            return '-1',[-1,]
296        else:
297            return '1',[1,]
298    elif SGData['SGLaue'] == '2/m': #monoclinic - default for 2D modulation vector
299        if '/' in SGData['SpGrp']:
300            return '2/m',[-1,1]
301        elif '2' in SGData['SpGrp']:
302            return '2',[-1,]
303        else:
304            return 'm',[1,]
305    elif SGData['SGLaue'] == 'mmm': #orthorhombic
306        if SGData['SpGrp'].count('2') == 3:
307            return '222',[-1,-1,-1]
308        elif SGData['SpGrp'].count('2') == 1:
309            if SGData['SGPolax'] == 'x':
310                return '2mm',[-1,1,1]
311            elif SGData['SGPolax'] == 'y':
312                return 'm2m',[1,-1,1]
313            elif SGData['SGPolax'] == 'z':
314                return 'mm2',[1,1,-1]
315        else:
316            return 'mmm',[1,1,1]
317    elif SGData['SGLaue'] == '4/m': #tetragonal
318        if '/' in SGData['SpGrp']:
319            return '4/m',[1,-1]
320        elif '-' in Flds[1]:
321            return '-4',[-1,]
322        else:
323            return '4',[1,]
324    elif SGData['SGLaue'] == '4/mmm':
325        if '/' in SGData['SpGrp']:
326            return '4/mmm',[1,-1,1,1]
327        elif '-' in Flds[1]:
328            if '2' in Flds[2]:
329                return '-42m',[-1,-1,1]
330            else:
331                return '-4m2',[-1,1,-1]             
332        elif '2' in Flds[2:]:
333            return '422',[1,-1,-1]
334        else:
335            return '4mm',[1,1,1]
336    elif SGData['SGLaue'] in ['3','3R']:  #trigonal/rhombohedral
337        if '-' in Flds[1]:
338            return '-3',[-1,]
339        else:
340            return '3',[1,]
341    elif SGData['SGLaue'] == '3mR' or 'R' in Flds[0]:
342        if '2' in Flds[2]:
343            return '32',[1,-1]
344        elif '-' in Flds[1]:
345            return '-3m',[-1,1]
346        else:
347            return '3m',[1,1]
348    elif SGData['SGLaue'] == '3m1':
349        if '2' in Flds[2]:
350            return '321',[1,-1,1]
351        elif '-' in Flds[1]:
352            return '-3m1',[-1,1,1]
353        else:
354            return '3m1',[1,1,1]
355    elif SGData['SGLaue'] == '31m':
356        if '2' in Flds[3]:
357            return '312',[1,1,-1]
358        elif '-' in Flds[1]:
359            return '-31m',[-1,1,1]
360        else:
361            return '31m',[1,1,1]
362    elif SGData['SGLaue'] == '6/m': #hexagonal
363        if '/' in SGData['SpGrp']:
364            return '6/m',[1,-1]
365        elif '-' in SGData['SpGrp']:
366            return '-6',[-1,]
367        else:
368            return '6',[1,]
369    elif SGData['SGLaue'] == '6/mmm':
370        if '/' in SGData['SpGrp']:
371            return '6/mmm',[1,-1,1,1]
372        elif '-' in Flds[1]:
373            if '2' in Flds[2]:
374                return '-62m',[-1,-1,1]
375            else:
376                return '-6m2',[-1,1,-1]                 
377        elif '2' in Flds[2:]:
378            return '622',[1,-1,-1]
379        else:
380            return '6mm',[1,1,1]   
381    elif SGData['SGLaue'] == 'm3':      #cubic - no (3+1) supersymmetry
382        if '2' in Flds[1]:
383            return '23',[]
384        else: 
385            return 'm3',[]
386    elif SGData['SGLaue'] == 'm3m':
387        if '4' in Flds[1]:
388            if '-' in Flds[1]:
389                return '-43m',[]
390            else:
391                return '432',[]
392        else:
393            return 'm3m',[]
394   
395def SGPrint(SGData,AddInv=False):
396    '''
397    Print the output of SpcGroup in a nicely formatted way. Used in SpaceGroup
398
399    :param SGData: from :func:`SpcGroup`
400    :returns:
401        SGText - list of strings with the space group details
402        SGTable - list of strings for each of the operations
403    '''
404    if SGData.get('SGFixed',False):       #inverses included in ops for cif fixed
405        Mult = len(SGData['SGCen'])*len(SGData['SGOps'])
406    else:
407        Mult = len(SGData['SGCen'])*len(SGData['SGOps'])*(int(SGData['SGInv'])+1)
408    SGText = []
409    SGText.append(' Space Group: '+SGData['SpGrp'])
410    if SGData.get('SGGray',False): SGText[-1] += " 1'"
411    CentStr = 'centrosymmetric'
412    if not SGData['SGInv']:
413        CentStr = 'non'+CentStr
414    if SGData['SGLatt'] in 'ABCIFR':
415        SGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
416    else:
417        SGText.append(' The lattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower()) 
418    SGText.append(' The Laue symmetry is '+SGData['SGLaue'])
419    if 'SGPtGrp' in SGData:         #patch
420        SGText.append(' The lattice point group is '+SGData['SGPtGrp'])
421    SGText.append(' Multiplicity of a general site is '+str(Mult))
422    if SGData['SGUniq'] in ['a','b','c']:
423        SGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
424    if SGData['SGInv']:
425        SGText.append(' The inversion center is located at 0,0,0')
426    if SGData['SGPolax']:
427        SGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
428    SGText.append(' ')
429    if len(SGData['SGCen']) == 1:
430        SGText.append(' The equivalent positions are:\n')
431    else:   
432        SGText.append(' The equivalent positions are:\n')
433        SGText.append(' ('+Latt2text(SGData['SGCen'])+')+\n')
434    SGTable = []
435    for i,Opr in enumerate(SGData['SGOps']):
436        SGTable.append('(%2d) %s'%(i+1,MT2text(Opr)))
437    if AddInv and SGData['SGInv']:
438        for i,Opr in enumerate(SGData['SGOps']):
439            IOpr = [-Opr[0],-Opr[1]]
440            SGTable.append('(%2d) %s'%(i+1,MT2text(IOpr)))       
441#    if SGData.get('SGGray',False) and not SGData.get('SGFixed',False):
442#        SGTable.append("     for 1'")
443#        for i,Opr in enumerate(SGData['SGOps']):
444#            SGTable.append('(%2d) %s'%(i+1,MT2text(Opr)))
445#        if AddInv and SGData['SGInv']:
446#            for i,Opr in enumerate(SGData['SGOps']):
447#                IOpr = [-Opr[0],-Opr[1]]
448#                SGTable.append('(%2d) %s'%(i+1,MT2text(IOpr)))       
449    return SGText,SGTable
450
451def AllOps(SGData):
452    '''
453    Returns a list of all operators for a space group, including those for
454    centering and a center of symmetry
455   
456    :param SGData: from :func:`SpcGroup`
457    :returns: (SGTextList,offsetList,symOpList,G2oprList) where
458
459      * SGTextList: a list of strings with formatted and normalized
460        symmetry operators.
461      * offsetList: a tuple of (dx,dy,dz) offsets that relate the GSAS-II
462        symmetry operation to the operator in SGTextList and symOpList.
463        these dx (etc.) values are added to the GSAS-II generated
464        positions to provide the positions that are generated
465        by the normalized symmetry operators.       
466      * symOpList: a list of tuples with the normalized symmetry
467        operations as (M,T) values
468        (see ``SGOps`` in the :ref:`Space Group object<SGData_table>`)
469      * G2oprList: The GSAS-II operations for each symmetry operation as
470        a tuple with (center,mult,opnum,opcode), where center is (0,0,0), (0.5,0,0),
471        (0.5,0.5,0.5),...; where mult is 1 or -1 for the center of symmetry
472        where opnum is the number for the symmetry operation, in ``SGOps``
473        (starting with 0) and opcode is mult*(100*icen+j+1).
474    '''
475    SGTextList = []
476    offsetList = []
477    symOpList = []
478    G2oprList = []
479    G2opcodes = []
480    onebar = (1,)
481    if SGData['SGInv']:
482        onebar += (-1,)
483    for icen,cen in enumerate(SGData['SGCen']):
484        for mult in onebar:
485            for j,(M,T) in enumerate(SGData['SGOps']):
486                offset = [0,0,0]
487                Tprime = (mult*T)+cen
488                for i in range(3):
489                    while Tprime[i] < 0:
490                        Tprime[i] += 1
491                        offset[i] += 1
492                    while Tprime[i] >= 1:
493                        Tprime[i] += -1
494                        offset[i] += -1
495                Opr = [mult*M,Tprime]
496                OPtxt = MT2text(Opr)
497                SGTextList.append(OPtxt.replace(' ',''))
498                offsetList.append(tuple(offset))
499                symOpList.append((mult*M,Tprime))
500                G2oprList.append((cen,mult,j))
501                G2opcodes.append(mult*(100*icen+j+1))
502    return SGTextList,offsetList,symOpList,G2oprList,G2opcodes
503
504def TextOps(text,table,reverse=False):
505    ''' Makes formatted operator list
506        :param text,table: arrays of text made by SGPrint
507        :param reverse: True for x+1/2 form; False for 1/2+x form
508        :returns: OpText: full list of symmetry operators; one operation per line
509        generally printed to console for use via cut/paste in other programs, but
510        could be used for direct input
511    '''
512    OpText = []
513    Inv = True
514    if 'noncentro' in text[1]:
515        Inv = False
516    Cent = [[0,0,0],]
517    if '0,0,0' in text[-1]:
518        Cent = np.array(eval(text[-1].split('+')[0].replace(';','),(')))
519    OpsM = []
520    OpsT = []
521    for item in table:
522        M,T = Text2MT(item.split(')')[1].replace(' ',''),CIF=True)
523        OpsM.append(M)
524        OpsT.append(T)
525    OpsM = np.array(OpsM)
526    OpsT = np.array(OpsT)
527    if Inv:
528        OpsM = np.concatenate((OpsM,-OpsM))
529        OpsT = np.concatenate((OpsT,-OpsT%1.))
530    for cent in Cent:
531        for iop,opM in enumerate(list(OpsM)):
532            txt = MT2text([opM,(OpsT[iop]+cent)%1.],reverse)
533            OpText.append(txt.replace(' ','').lower())
534    return OpText
535
536def TextGen(SGData,reverse=False):
537    GenSym,GenFlg,BNSsym = GetGenSym(SGData)
538    SGData['GenSym'] = GenSym
539    SGData['GenFlg'] = GenFlg
540    text,table = SGPrint(SGData)
541    GenText = []
542    OprNames = GetOprNames(SGData)
543    OpText = TextOps(text,table,reverse)
544    for name in SGData['GenSym']:
545        gid = OprNames.index(name.replace(' ',''))
546        GenText.append(OpText[gid])
547    if len(SGData['SGCen']) > 1:
548        GenText.append(OpText[-1])
549    return GenText
550
551def GetOprNames(SGData):
552    OprNames = [GetOprPtrName(str(irtx)) for irtx in PackRot(SGData['SGOps'])]
553    if SGData['SGInv']:
554        OprNames += [GetOprPtrName(str(-irtx)) for irtx in PackRot(SGData['SGOps'])]
555    return OprNames
556   
557def MT2text(Opr,reverse=False):
558    "From space group matrix/translation operator returns text version"
559    XYZ = ('-Z','-Y','-X','X-Y','ERR','Y-X','X','Y','Z')
560    TRA = ('   ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
561    Fld = ''
562    M,T = Opr
563    for j in range(3):
564        IJ = int(round(2*M[j][0]+3*M[j][1]+4*M[j][2]+4))%12
565        IK = int(round(T[j]*12))%12
566        if IK:
567            if IJ < 3:
568                if reverse:
569                    Fld += (XYZ[IJ]+'+'+TRA[IK]).rjust(5)
570                else:
571                    Fld += (TRA[IK]+XYZ[IJ]).rjust(5)
572            else:
573                if reverse:
574                    Fld += (XYZ[IJ]+'+'+TRA[IK]).rjust(5)
575                else:
576                    Fld += (TRA[IK]+'+'+XYZ[IJ]).rjust(5)
577        else:
578            Fld += XYZ[IJ].rjust(5)
579        if j != 2: Fld += ', '
580    return Fld
581   
582def Latt2text(Cen):
583    "From lattice centering vectors returns ';' delimited cell centering vectors"
584    lattTxt = ''
585    fracList = ['1/2','1/3','2/3','1/4','3/4','1/5','2/5','3/5','4/5','1/6','5/6',
586        '1/7','2/7','3/7','4/7','5/7','6/7','1/8','3/8','5/8','7/8','1/9','2/9','4/9','5/9','7/9','8/9']
587    mulList = [2,3,3,4,4,5,5,5,5,6,6,7,7,7,7,7,7,8,8,8,8,9,9,9,9,9,9]
588    prodList = [1.,1.,2.,1.,3.,1.,2.,3.,4.,1.,5.,1.,2.,3.,4.,5.,6.,1.,3.,5.,7.,1.,2.,4.,5.,7.,8.]
589    nCen = len(Cen)
590    for i,cen in enumerate(Cen):
591        txt = ''
592        for icen in cen:
593            if icen == 1:
594                txt += '1,'
595                continue
596            if not icen:
597                txt += '0,'
598                continue
599            if icen < 0:
600                txt += '-'
601                icen *= -1
602            for mul,prod,frac in zip(mulList,prodList,fracList):
603                if abs(icen*mul-prod) < 1.e-5:
604                    txt += frac+','
605                    break
606        lattTxt += txt[:-1]+'; '
607        if i and not i%8 and i < nCen-1:    #not for the last cen!
608            lattTxt += '\n     '
609    return lattTxt[:-2]
610
611def SpaceGroup(SGSymbol):
612    '''
613    Print the output of SpcGroup in a nicely formatted way.
614
615    :param SGSymbol: space group symbol (string) with spaces between axial fields
616    :returns: nothing
617    '''
618    E,A = SpcGroup(SGSymbol)
619    if E > 0:
620        print (SGErrors(E))
621        return
622    for l in SGPrint(A):
623        print (l)
624################################################################################
625#### Magnetic space group stuff
626################################################################################
627
628def GetGenSym(SGData):
629    '''
630    Get the space group generator symbols
631    :param SGData: from :func:`SpcGroup`
632    LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
633    LattSym = ('P','A','B','C','I','F','R')
634   
635    '''
636    OprNames = [GetOprPtrName(str(irtx)) for irtx in PackRot(SGData['SGOps'])]
637    if SGData['SGInv']:
638        OprNames += [GetOprPtrName(str(-irtx)) for irtx in PackRot(SGData['SGOps'])]
639    Nsyms = len(SGData['SGOps'])
640    if SGData['SGInv'] and not SGData['SGFixed']: Nsyms *= 2
641    UsymOp = ['1',]
642    OprFlg = [0,] 
643    if Nsyms == 2:                    #Centric triclinic or acentric monoclinic
644        UsymOp.append(OprNames[1])
645        OprFlg.append(SGData['SGGen'][1])
646    elif Nsyms == 4:                    #Point symmetry 2/m, 222, 22m, or 4
647        if '4z' in OprNames[1]:          #Point symmetry 4 or -4
648            UsymOp.append(OprNames[1])
649            OprFlg.append(SGData['SGGen'][1])
650        elif not SGData['SGInv']:       #Acentric Orthorhombic
651            if 'm' in OprNames[1:4]:    #22m, 2m2 or m22
652                if '2' in OprNames[1]:      #Acentric orthorhombic, 2mm
653                    UsymOp.append(OprNames[2])
654                    OprFlg.append(SGData['SGGen'][2])
655                    UsymOp.append(OprNames[3])
656                    OprFlg.append(SGData['SGGen'][3])
657                elif '2' in OprNames[2]:    #Acentric orthorhombic, m2m
658                    UsymOp.append(OprNames[1])
659                    OprFlg.append(SGData['SGGen'][1])
660                    UsymOp.append(OprNames[3])
661                    OprFlg.append(SGData['SGGen'][3])
662                else:                       #Acentric orthorhombic, mm2
663                    UsymOp.append(OprNames[1])
664                    OprFlg.append(SGData['SGGen'][1])
665                    UsymOp.append(OprNames[2])
666                    OprFlg.append(SGData['SGGen'][2])
667            else:                           #Acentric orthorhombic, 222
668                SGData['SGGen'][1:] = [4,2,1]
669                UsymOp.append(OprNames[1])
670                OprFlg.append(SGData['SGGen'][1])
671                UsymOp.append(OprNames[2])
672                OprFlg.append(SGData['SGGen'][2])
673                UsymOp.append(OprNames[3])
674                OprFlg.append(SGData['SGGen'][3])
675        else:                               #Centric Monoclinic
676            UsymOp.append(OprNames[1])
677            OprFlg.append(SGData['SGGen'][1])
678            UsymOp.append(OprNames[3])
679            OprFlg.append(SGData['SGGen'][3])
680    elif Nsyms == 6:                    #Point symmetry 32, 3m or 6
681            if '6' in OprNames[1]:      #Hexagonal 6/m Laue symmetry
682                UsymOp.append(OprNames[1])
683                OprFlg.append(SGData['SGGen'][1])
684            else:                       #Trigonal
685                UsymOp.append(OprNames[4])
686                OprFlg.append(SGData['SGGen'][3])
687                if '2110' in OprNames[1]: UsymOp[-1] = ' 2100 '
688    elif Nsyms == 8:                    #Point symmetry mmm, 4/m, or 422, etc
689        if '4' in OprNames[1]:           #Tetragonal
690            if SGData['SGInv']:         #4/m
691                UsymOp.append(OprNames[1])
692                OprFlg.append(SGData['SGGen'][1])
693                UsymOp.append(OprNames[6])
694                OprFlg.append(SGData['SGGen'][6])
695            else:
696                if 'x' in OprNames[4]:      #4mm type group
697                    UsymOp.append(OprNames[4])
698                    OprFlg.append(6)
699                    UsymOp.append(OprNames[7])
700                    OprFlg.append(8)
701                else:                       #-42m, -4m2, and 422 type groups
702                    UsymOp.append(OprNames[5])
703                    OprFlg.append(8)
704                    UsymOp.append(OprNames[6])
705                    OprFlg.append(19)
706        else:                               #Orthorhombic, mmm
707            UsymOp.append(OprNames[1])
708            OprFlg.append(SGData['SGGen'][1])
709            UsymOp.append(OprNames[2])
710            OprFlg.append(SGData['SGGen'][2])
711            UsymOp.append(OprNames[7])
712            OprFlg.append(SGData['SGGen'][7])
713    elif Nsyms == 12 and '3' in OprNames[1] and SGData['SGSys'] != 'cubic':        #Trigonal
714        UsymOp.append(OprNames[3])
715        OprFlg.append(SGData['SGGen'][3])
716        UsymOp.append(OprNames[9])
717        OprFlg.append(SGData['SGGen'][9])
718    elif Nsyms == 12 and '6' in OprNames[1]:        #Hexagonal
719        if 'mz' in OprNames[9]:                     #6/m
720            UsymOp.append(OprNames[1])
721            OprFlg.append(SGData['SGGen'][1])
722            UsymOp.append(OprNames[6])
723            OprFlg.append(SGData['SGGen'][6])
724        else:                                       #6mm, -62m, -6m2 or 622
725            UsymOp.append(OprNames[6])
726            OprFlg.append(18)
727            if 'm' in UsymOp[-1]: OprFlg[-1] = 20
728            UsymOp.append(OprNames[7])
729            OprFlg.append(24)
730    elif Nsyms in [16,24]:
731        if '3' in OprNames[1]:
732            UsymOp.append('')
733            OprFlg.append(SGData['SGGen'][3])
734            for i in range(Nsyms):
735                if 'mx' in OprNames[i]:
736                    UsymOp[-1] = OprNames[i]
737                elif 'm11' in OprNames[i]:
738                    UsymOp[-1] = OprNames[i]
739                elif '211' in OprNames[i]:
740                    UsymOp[-1] = OprNames[i]
741                    OprFlg[-1] = 24
742        else:                                     #4/mmm or 6/mmm
743            UsymOp.append('  mz  ')
744            OprFlg.append(1)
745            if '4' in OprNames[1]:                  #4/mmm
746                UsymOp.append('  mx  ')
747                OprFlg.append(20)
748                UsymOp.append(' m110 ')
749                OprFlg.append(24)
750            else:                                   #6/mmm
751                UsymOp.append(' m110 ')
752                OprFlg.append(4)
753                UsymOp.append(' m+-0 ')
754                OprFlg.append(8)
755    else:                                           #System is cubic
756        if Nsyms == 48:
757            UsymOp.append('  mx  ')
758            OprFlg.append(4)
759            UsymOp.append(' m110 ')
760            OprFlg.append(24)
761           
762    if 'P' in SGData['SGLatt']:
763        if SGData['SGSys'] == 'triclinic':
764            BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5]}           
765        elif SGData['SGSys'] == 'monoclinic':
766            BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5]}
767            if SGData['SGUniq'] == 'a':
768                BNSsym.update({'P_B':[.5,0,.5],'P_C':[.5,.5,0]})
769            elif SGData['SGUniq'] == 'b':
770                BNSsym.update({'P_A':[.5,.5,0],'P_C':[0,.5,.5]})
771            elif SGData['SGUniq'] == 'c':
772                BNSsym.update({'P_A':[0,.5,.5],'P_B':[.5,0,.5]})
773        elif SGData['SGSys'] == 'orthorhombic':
774            BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5],
775                'P_A':[0,.5,.5],'P_B':[.5,0,.5],'P_C':[.5,.5,0],'P_I':[.5,.5,.5]}
776        elif SGData['SGSys'] == 'tetragonal':
777            BNSsym = {'P_c':[0,0,.5],'P_C':[.5,.5,0],'P_I':[.5,.5,.5]}           
778        elif SGData['SGSys'] in ['trigonal','hexagonal']:
779            BNSsym = {'P_c':[0,0,.5]}           
780        elif SGData['SGSys'] == 'cubic':
781            BNSsym = {'P_I':[.5,.5,.5]}           
782           
783    elif 'A' in SGData['SGLatt']:
784        if SGData['SGSys'] == 'monoclinic':
785            BNSsym = {}
786            if SGData['SGUniq'] == 'b':
787                BNSsym.update({'A_a':[.5,0,0],'A_c':[0,0,.5]})
788            elif SGData['SGUniq'] == 'c':
789                BNSsym.update({'A_a':[.5,0,0],'A_b':[0,.5,0]})
790        elif SGData['SGSys'] == 'orthorhombic':
791            BNSsym = {'A_a':[.5,0,0],'A_b':[0,.5,0],'A_c':[0,0,.5],
792               'A_B':[.5,0,.5],'A_C':[.5,.5,0]}   
793        elif SGData['SGSys'] == 'triclinic':
794            BNSsym = {'A_a':[.5,0,0],'A_b':[0,.5,0],'A_c':[0,0,.5]}   
795           
796    elif 'B' in SGData['SGLatt']:
797        if SGData['SGSys'] == 'monoclinic':
798            BNSsym = {}
799            if SGData['SGUniq'] == 'a':
800                BNSsym.update({'B_b':[0,.5,0],'B_c':[0,0,.5]})
801            elif SGData['SGUniq'] == 'c':
802                BNSsym.update({'B_a':[.5,0,0],'B_b':[0,.5,0]})
803        elif SGData['SGSys'] == 'orthorhombic':
804            BNSsym = {'B_a':[.5,0,0],'B_b':[0,.5,0],'B_c':[0,0,.5],
805                'B_A':[0,.5,.5],'B_C':[.5,.5,0]}     
806        elif SGData['SGSys'] == 'triclinic':
807            BNSsym = {'B_a':[.5,0,0],'B_b':[0,.5,0],'B_c':[0,0,.5]}     
808           
809    elif 'C' in SGData['SGLatt']:
810        if SGData['SGSys'] == 'monoclinic':
811            BNSsym = {}
812            if SGData['SGUniq'] == 'a':
813                BNSsym.update({'C_b':[0,.5,.0],'C_c':[0,0,.5]})
814            elif SGData['SGUniq'] == 'b':
815                BNSsym.update({'C_a':[.5,0,0],'C_c':[0,0,.5]})
816        elif SGData['SGSys'] == 'orthorhombic':
817            BNSsym = {'C_a':[.5,0,0],'C_b':[0,.5,0],'C_c':[0,0,.5],
818                'C_A':[0,.5,.5],'C_B':[.5,0,.5]}     
819        elif SGData['SGSys'] == 'triclinic':
820            BNSsym = {'C_a':[.5,0,0],'C_b':[0,.5,0],'C_c':[0,0,.5]}     
821           
822    elif 'I' in SGData['SGLatt']:
823        if SGData['SGSys'] in ['monoclinic','orthorhombic','triclinic']:
824            BNSsym = {'I_a':[.5,0,0],'I_b':[0,.5,0],'I_c':[0,0,.5]}
825        elif SGData['SGSys'] == 'tetragonal':
826            BNSsym = {'I_c':[0,0,.5]}
827        elif SGData['SGSys'] == 'cubic':
828            BNSsym = {} 
829           
830    elif 'F' in SGData['SGLatt']:
831        if SGData['SGSys'] in ['monoclinic','orthorhombic','cubic','triclinic']:
832            BNSsym = {'F_S':[.5,.5,.5]}
833           
834    elif 'R' in SGData['SGLatt']:
835        BNSsym = {'R_I':[0,0,.5]}
836    return UsymOp,OprFlg,BNSsym
837
838def ApplyBNSlatt(SGData,BNSlatt):
839    Tmat = np.eye(3)
840    BNS = BNSlatt[0]
841    A = np.array(BNSlatt[1])
842    Laue = SGData['SGLaue']
843    SGCen = SGData['SGCen']
844    if '_a' in BNS:
845        Tmat[0,0] = 2.0
846    elif '_b' in BNS:
847        Tmat[1,1] = 2.0
848    elif '_c' in BNS:
849        Tmat[2,2] = 2.0
850    elif '_A' in BNS:
851        Tmat[0,0] = 2.0
852    elif '_B' in BNS:
853        Tmat[1,1] = 2.0
854    elif '_C' in BNS:
855        Tmat[2,2] = 2.0
856    elif '_I' in BNS:
857        Tmat *= 2.0
858        if 'R' in Laue:
859            SGData['SGSpin'][-1] = -1
860        else:
861            SGData['SGSpin'].append(-1)
862    elif '_S' in BNS:
863        SGData['SGSpin'][-1] = -1
864        SGData['SGSpin'] += [-1,-1,-1,]
865        Tmat *= 2.0
866    else:
867        return Tmat
868    SGData['SGSpin'].append(-1)
869    C = SGCen+A
870    SGData['SGCen'] = np.vstack((SGCen,C))%1.
871    return Tmat
872       
873def CheckSpin(isym,SGData):
874    ''' Check for exceptions in spin rules
875    '''
876    if SGData['SGPtGrp'] in ['222','mm2','2mm','m2m']:      #only 2/3 can be red; not 1/3 or 3/3
877        if SGData['SGSpin'][1]*SGData['SGSpin'][2]*SGData['SGSpin'][3] < 0:
878            SGData['SGSpin'][(isym+1)%3+1] *= -1
879        if SGData['SpGrp'][0] == 'F' and isym > 2:
880            SGData['SGSpin'][(isym+1)%3+3] == 1
881    elif SGData['SGPtGrp'] == 'mmm':
882        if SGData['SpGrp'][0] == 'F' and isym > 2:
883            SGData['SGSpin'][(isym+1)%3+3] == 1
884
885def MagSGSym(SGData):       #needs to use SGPtGrp not SGLaue!
886    SGLaue = SGData['SGLaue']
887    if '1' not in SGData['GenSym']:        #patch for old gpx files
888        SGData['GenSym'] = ['1',]+SGData['GenSym']
889        SGData['SGSpin'] = [1,]+list(SGData['SGSpin'])
890    if len(SGData['SGSpin'])<len(SGData['GenSym']):
891        SGData['SGSpin'] = [1,]+list(SGData['SGSpin'])      #end patch
892    GenSym = SGData['GenSym'][1:]       #skip identity
893    SpnFlp = SGData['SGSpin']
894#    print('SpnFlp',SpnFlp)
895    SGPtGrp = SGData['SGPtGrp']
896    if len(SpnFlp) == 1:
897        SGData['MagPtGp'] = SGPtGrp
898        return SGData['SpGrp']
899    magSym = SGData['SpGrp'].split()
900    if SGLaue in ['-1',]:
901        SGData['MagPtGp'] = SGPtGrp
902        if SpnFlp[1] == -1:
903            magSym[1] += "'"
904            SGData['MagPtGp'] += "'"
905    elif SGLaue in ['2/m','4/m','6/m']: #all ok
906        Uniq = {'a':1,'b':2,'c':3,'':1}
907        id = [0,1]
908        if len(magSym) > 2:
909            id = [0,Uniq[SGData['SGUniq']]]
910        sym = magSym[id[1]].split('/')
911        Ptsym = SGLaue.split('/')
912        if len(GenSym) == 3:
913            for i in [0,1,2]:
914                if SpnFlp[i+1] < 0:
915                    sym[i] += "'"
916                    Ptsym[i] += "'"
917        else:
918            for i in range(len(GenSym)):
919                if SpnFlp[i+1] < 0:                     
920                    sym[i] += "'"
921                    Ptsym[i] += "'"
922        SGData['MagPtGp'] = '/'.join(Ptsym)
923        magSym[id[1]] = '/'.join(sym)
924    elif SGPtGrp in ['mmm','mm2','m2m','2mm','222']:
925        SGData['MagPtGp'] = ''
926        for i in [0,1,2]:
927            SGData['MagPtGp'] += SGPtGrp[i]
928            if SpnFlp[i+1] < 0:
929                magSym[i+1] += "'"
930                SGData['MagPtGp'] += "'"
931    elif SGLaue == '6/mmm': #ok
932        magPtGp = list(SGPtGrp)
933        if len(GenSym) == 2:
934            for i in [0,1]:
935                if SpnFlp[i+1] < 0:
936                    magSym[i+2] += "'"
937                    magPtGp[i+1] += "'"
938            if SpnFlp[1]*SpnFlp[2] < 0:
939                magSym[1] += "'"
940                magPtGp[0] += "'"
941        else:
942            sym = magSym[1].split('/')
943            Ptsym = ['6','m']
944            magPtGp = ['','m','m']
945            for i in [0,1,2]:
946                if SpnFlp[i+1] < 0:
947                    if i:
948                        magSym[i+1] += "'"
949                        magPtGp[i] += "'"
950                    else:
951                        sym[1] += "'"
952                        Ptsym[0] += "'"
953            if SpnFlp[2]*SpnFlp[3] < 0:
954                sym[0] += "'"                   
955                Ptsym[0] += "'"                   
956            magSym[1] = '/'.join(sym)
957            magPtGp[0] = '/'.join(Ptsym)
958        SGData['MagPtGp'] = ''.join(magPtGp)
959    elif SGLaue == '4/mmm':
960        magPtGp = list(SGPtGrp)
961        if len(GenSym) == 2:
962            for i in [0,1]:
963                if SpnFlp[i+1] < 0:
964                    magSym[i+2] += "'"
965                    magPtGp[i+1] += "'"
966            if SpnFlp[1]*SpnFlp[2] < 0:
967                magSym[1] += "'"
968                magPtGp[0] += "'"
969        else:
970            if '/' in magSym[1]:    #P 4/m m m, etc.
971                sym = magSym[1].split('/')
972                Ptsym = ['4','m']
973                magPtGp = ['','m','m']
974                for i in [0,1,2]:
975                    if SpnFlp[i+1] < 0:
976                        if i:
977                            magSym[i+1] += "'"
978                            magPtGp[i] += "'"
979                        else:
980                            sym[1] += "'"
981                            Ptsym[1] += "'"
982                if SpnFlp[2]*SpnFlp[3] < 0:
983                    sym[0] += "'"                   
984                    Ptsym[0] += "'"                   
985                magSym[1] = '/'.join(sym)
986                magPtGp[0] = '/'.join(Ptsym)
987            else:
988                for i in [0,1]:
989                    if SpnFlp[i+1] < 0:
990                        magSym[i+2] += "'"
991                if SpnFlp[1]*SpnFlp[2] < 0:
992                    magSym[1] += "'"
993        SGData['MagPtGp'] = ''.join(magPtGp)
994    elif SGLaue in ['3','3m1','31m']:   #ok
995        Ptsym = list(SGLaue)
996        if len(GenSym) == 1:    #all ok
997            id = 2
998            if (len(magSym) == 4) and (magSym[2] == '1'):
999                id = 3
1000            if '3' in GenSym[0]:
1001                id = 1
1002            magSym[id].strip("'")
1003            if SpnFlp[1] < 0:
1004                magSym[id] += "'"
1005                Ptsym[id-1] += "'"
1006        elif len(GenSym) == 2:
1007            if 'R' in GenSym[1]:
1008                magSym[-1].strip("'")
1009                if SpnFlp[1] < 0:
1010                    magSym[-1] += "'"
1011                    Ptsym[-1] += "'"
1012            else:
1013                i,j = [1,2]
1014                if magSym[2] == '1':
1015                    i,j = [1,3]
1016                magSym[i].strip("'")
1017                Ptsym[i-1].strip("'")
1018                magSym[j].strip("'")
1019                Ptsym[j-1].strip("'")
1020                if SpnFlp[1:3] == [1,-1]:
1021                    magSym[i] += "'"
1022                    Ptsym[i-1] += "'"
1023                elif SpnFlp[1:3] == [-1,-1]:
1024                    magSym[j] += "'"
1025                    Ptsym[j-1] += "'"
1026                elif SpnFlp[1:3] == [-1,1]:
1027                    magSym[i] += "'"
1028                    Ptsym[i-1] += "'"
1029                    magSym[j] += "'"
1030                    Ptsym[j-1] += "'"
1031        elif len(GenSym):
1032            if 'c' not in magSym[2]:
1033                i,j = [1,2]
1034                magSym[i].strip("'")
1035                Ptsym[i-1].strip("'")
1036                magSym[j].strip("'")
1037                Ptsym[j-1].strip("'")
1038                if SpnFlp[1:3] == [1,-1]:
1039                    magSym[i] += "'"
1040                    Ptsym[i-1] += "'"
1041                elif SpnFlp[1:3] == [-1,-1]:
1042                    magSym[j] += "'"
1043                    Ptsym[j-1] += "'"
1044                elif SpnFlp[2] == [-1,1]:
1045                    magSym[i] += "'"
1046                    Ptsym[i-1] += "'"
1047                    magSym[j] += "'"
1048                    Ptsym[j-1] += "'"
1049        SGData['MagPtGp'] = ''.join(Ptsym)
1050    elif SGData['SGPtGrp'] == '23' and len(magSym):
1051        SGData['MagPtGp'] = '23'
1052    elif SGData['SGPtGrp'] == 'm3':
1053        SGData['MagPtGp'] = "m3"
1054        if SpnFlp[1] < 0:
1055            magSym[1] += "'"
1056            magSym[2] += "'"
1057            SGData['MagPtGp'] = "m'3'"
1058        if SpnFlp[2] < 0:
1059            if not 'm' in magSym[1]:    #only Ia3
1060                magSym[1].strip("'")
1061                SGData['MagPtGp'] = "m3'"
1062    elif SGData['SGPtGrp'] in ['432','-43m']:
1063        Ptsym = SGData['SGPtGrp'].split('3')
1064        if SpnFlp[1] < 0:
1065            magSym[1] += "'"
1066            Ptsym[0] += "'"
1067            magSym[3] += "'"
1068            Ptsym[1] += "'"
1069        SGData['MagPtGp'] = '3'.join(Ptsym)
1070    elif SGData['SGPtGrp'] == 'm3m':
1071        Ptsym = ['m','3','m']
1072        if SpnFlp[1:3] == [-1,1]:
1073            magSym[1] += "'"
1074            Ptsym[0] += "'"
1075            magSym[2] += "'"
1076            Ptsym[1] += "'"
1077        elif SpnFlp[1:3] == [1,-1]:
1078            magSym[3] += "'"
1079            Ptsym[2] += "'"
1080        elif SpnFlp[1:3] == [-1,-1]:
1081            magSym[1] += "'"
1082            Ptsym[0] += "'"
1083            magSym[2] += "'"
1084            Ptsym[1] += "'"
1085            magSym[3] += "'"
1086            Ptsym[2] += "'"
1087        SGData['MagPtGp'] = ''.join(Ptsym)
1088#    print SpnFlp
1089    magSym[0] = SGData.get('BNSlattsym',[SGData['SGLatt'],[0,0,0]])[0]
1090    return ' '.join(magSym)
1091
1092def fixMono(SpGrp):
1093    'fixes b-unique monoclinics in e.g. P 1 2/1c 1 --> P 21/c '
1094    Flds = SpGrp.split()
1095    if len(Flds) == 4:
1096        if Flds[2] != '1':
1097            return '%s %s'%(Flds[0],Flds[2])
1098        else:
1099            return None
1100    else:
1101        return SpGrp
1102
1103def Trans2Text(Trans):
1104    "from transformation matrix to text"
1105    cells = ['a','b','c']
1106    Text = ''
1107    for row in Trans:
1108        Fld = ''
1109        for i in [0,1,2]:
1110            if row[i]:
1111                if Fld and row[i] > 0.:
1112                    Fld += '+'
1113                Fld += '%3.1f'%(row[i])+cells[i]
1114        Text += Fld
1115        Text += ','
1116        Text = Text.replace('1.0','').replace('.0','').replace('0.5','1/2')
1117    return Text[:-1]
1118
1119def getlattSym(Trans):
1120    Fives = {'ababc':'abc','bcbca':'cba','acacb':'acb','cabab':'cab','abcab':'acb'}
1121    transText = Trans2Text(Trans)
1122    lattSym = ''
1123    for fld in transText.split(','):
1124        if 'a' in fld: lattSym += 'a'
1125        if 'b' in fld: lattSym += 'b'
1126        if 'c' in fld: lattSym += 'c'
1127    if len(lattSym) == 5:
1128        lattSym = Fives[lattSym]
1129    return lattSym
1130
1131
1132def Text2MT(mcifOpr,CIF=True):
1133    "From space group cif text returns matrix/translation"
1134    XYZ = {'x':[1,0,0],'+x':[1,0,0],'-x':[-1,0,0],'y':[0,1,0],'+y':[0,1,0],'-y':[0,-1,0],
1135           'z':[0,0,1],'+z':[0,0,1],'-z':[0,0,-1],'x-y':[1,-1,0],'-x+y':[-1,1,0],'y-x':[-1,1,0],
1136           '+x-y':[1,-1,0],'+y-x':[-1,1,0]}
1137    ops = mcifOpr.split(",")
1138    M = []
1139    T = []
1140    for op in ops[:3]:
1141        ip = len(op)
1142        if '/' in op:
1143            try:    #mcif format
1144                nP = op.count('+')
1145                opMT = op.split('+')
1146                T.append(eval(opMT[nP]))
1147                if nP == 2:
1148                    opMT[0] = '+'.join(opMT[0:2])
1149            except NameError:   #normal cif format
1150                ip = op.index('/')
1151                T.append(eval(op[:ip+2]))
1152                opMT = [op[ip+2:],'']
1153        else:
1154            opMT = [op,'']
1155            T.append(0.)
1156        M.append(XYZ[opMT[0].lower()])
1157    return np.array(M),np.array(T)
1158           
1159def MagText2MTS(mcifOpr,CIF=True):
1160    "From magnetic space group cif text returns matrix/translation + spin flip"
1161    XYZ = {'x':[1,0,0],'+x':[1,0,0],'-x':[-1,0,0],'y':[0,1,0],'+y':[0,1,0],'-y':[0,-1,0],
1162           'z':[0,0,1],'+z':[0,0,1],'-z':[0,0,-1],'x-y':[1,-1,0],'-x+y':[-1,1,0],'y-x':[-1,1,0],
1163           '+x-y':[1,-1,0],'+y-x':[-1,1,0]}
1164    ops = mcifOpr.split(",")
1165    M = []
1166    T = []
1167    for op in ops[:3]:
1168        ip = len(op)
1169        if '/' in op:
1170            try:    #mcif format
1171                nP = op.count('+')
1172                opMT = op.split('+')
1173                T.append(eval(opMT[nP]))
1174                if nP == 2:
1175                    opMT[0] = '+'.join(opMT[0:2])
1176            except NameError:   #normal cif format
1177                ip = op.index('/')
1178                T.append(eval(op[:ip+2]))
1179                opMT = [op[ip+2:],'']
1180        else:
1181            opMT = [op,'']
1182            T.append(0.)
1183        M.append(XYZ[opMT[0].lower()])
1184    spnflp = 1
1185    if '-1' in ops[3]:
1186        spnflp = -1
1187    return np.array(M),np.array(T),spnflp
1188           
1189def MagSSText2MTS(mcifOpr):
1190    "From magnetic super space group cif text returns matrix/translation + spin flip"
1191    XYZ = {'x1':[1,0,0,0],'-x1':[-1,0,0,0],
1192           'x2':[0,1,0,0],'-x2':[0,-1,0,0],
1193           'x3':[0,0,1,0],'-x3':[0,0,-1,0],
1194           'x4':[0,0,0,1],'-x4':[0,0,0,-1],
1195           'x1-x2':[1,-1,0,0],'-x1+x2':[-1,1,0,0],
1196           'x1-x4':[1,0,0,-1],'-x1+x4':[-1,0,0,1],
1197           'x2-x4':[0,1,0,-1],'-x2+x4':[0,-1,0,1],
1198           '-x1-x2+x4':[-1,-1,0,1],'x1+x2-x4':[1,1,0,-1]}
1199    ops = mcifOpr.split(",")
1200    M = []
1201    T = []
1202    for op in ops[:4]:
1203        ip = len(op)
1204        if '/' in op:
1205            ip = op.index('/')-2
1206            T.append(eval(op[ip:]))
1207        else:
1208            T.append(0.)
1209        M.append(XYZ[op[:ip]])
1210    spnflp = 1
1211    if '-1' in ops[4]:
1212        spnflp = -1
1213    return np.array(M),np.array(T),spnflp
1214
1215def GetSGSpin(SGData,MSgSym):
1216    'get spin generators from magnetic space group symbol'
1217    SpGrp = SGData['SpGrp']
1218    mSgSym = MSgSym+' '
1219    Flds = SpGrp.split()
1220    iB = 0
1221    Spn = [1,]          #for identity generator
1222    if len(Flds) == 2:  #-1,  2/m, 4/m & 6/m; 1 or 2 generators
1223        fld = Flds[1]
1224        iF = mSgSym[iB:].index(fld[0])+iB
1225        jF = mSgSym[iF:].index(fld[-1])+iF
1226        if '/' in mSgSym[iF:jF]:
1227            if "'" in mSgSym[iF:jF]:
1228                Spn.append(-1)
1229            else:
1230                Spn.append(1)
1231        if "'" == mSgSym[jF+1]:
1232            Spn.append(-1)
1233        else:
1234            Spn.append(1)
1235    elif len(Flds) == 3:    # 3m & m3; 1 or 2 generator
1236        if SGData['SGPtGrp'] == '-3m':
1237            if not mSgSym.count("'"):
1238                Spn += [1,1,]
1239            elif mSgSym.count("'") == 2:
1240                Spn += [-1,1,]
1241            elif "3'" in mSgSym:
1242                Spn += [1,-1,]
1243            else:
1244                Spn += [-1,-1,]
1245        else:
1246            if "'" in mSgSym:   #could be 1 or 2 '; doesn't matter.
1247                Spn.append(-1)
1248            else:
1249                Spn.append(1)
1250    else:                   #the rest; 3 generators. NB:  any ' before / in 1st field ignored
1251        for fld in Flds[1:]:
1252            iF = mSgSym[iB:].index(fld[0])+iB
1253            jF = mSgSym[iF:].index(fld[-1])+iF
1254            if "'" == mSgSym[jF+1]:
1255                Spn.append(-1)
1256                iB = jF+2
1257            else:
1258                Spn.append(1)
1259                iB = jF+1
1260    return Spn
1261           
1262def GenMagOps(SGData):
1263    FlpSpn = SGData['SGSpin']
1264    Nsym = len(SGData['SGOps'])
1265    Ncv = len(SGData['SGCen'])
1266    sgOp = [M for M,T in SGData['SGOps']]
1267    oprName = [GetOprPtrName(str(irtx)) for irtx in PackRot(SGData['SGOps'])]
1268    if SGData['SGInv'] and not SGData['SGFixed'] and not SGData['SGGray']:
1269        Nsym *= 2
1270        sgOp += [-M for M,T in SGData['SGOps']]
1271        oprName += [GetOprPtrName(str(-irtx)) for irtx in PackRot(SGData['SGOps'])]
1272    Nsyms = 0
1273    sgOps = []
1274    OprNames = []
1275    for incv in range(Ncv):
1276        Nsyms += Nsym
1277        sgOps += sgOp
1278        OprNames += oprName
1279    if SGData['SGFixed']:
1280        SpnFlp = SGData['SpnFlp']
1281    else:
1282        SpnFlp = np.ones(Nsym,dtype=np.int)
1283        GenFlg = SGData.get('GenFlg',[0])
1284        Ngen = len(SGData['SGGen'])
1285    #    print ('GenFlg:',SGData['GenFlg'])
1286    #    print ('GenSym:',SGData['GenSym'])
1287        Nfl = len(GenFlg)
1288        for ieqv in range(Nsym):
1289            for iunq in range(Nfl):
1290                if SGData['SGGen'][ieqv%Ngen] & GenFlg[iunq]:
1291                    SpnFlp[ieqv] *= FlpSpn[iunq]
1292    #        print ('\nMagSpGrp:',SGData['MagSpGrp'],Ncv)
1293    #        print ('FlpSpn:',Nfl,FlpSpn)
1294        for incv in range(Ncv):
1295            if incv:
1296                try:
1297                    SpnFlp = np.concatenate((SpnFlp,SpnFlp[:Nsym]*FlpSpn[Nfl+incv-1]))
1298                except IndexError:
1299                    FlpSpn = [1,]+FlpSpn
1300                    SpnFlp = np.concatenate((SpnFlp,SpnFlp[:Nsym]*FlpSpn[Nfl+incv-1]))                   
1301    detM = [nl.det(M) for M in sgOp]
1302    MagMom = SpnFlp*np.array(Ncv*detM)      #duplicate for no. centerings
1303    SGData['MagMom'] = MagMom
1304#        print ('SgOps:',OprNames)
1305#        print ('SGGen:',SGData['SGGen'])
1306#        print ('SpnFlp:',SpnFlp)
1307#        print ('MagMom:',MagMom)
1308    return OprNames,SpnFlp
1309   
1310def GetOpNum(Opr,SGData):
1311    Nops = len(SGData['SGOps'])
1312    opNum = abs(Opr)%100
1313    cent = abs(Opr)//100
1314    if Opr < 0 and not SGData['SGFixed']:
1315        opNum += Nops
1316    if SGData['SGInv'] and not SGData['SGFixed']:
1317        Nops *= 2
1318    opNum += cent*Nops
1319    return opNum
1320       
1321################################################################################
1322#### Superspace group codes
1323################################################################################
1324       
1325def SSpcGroup(SGData,SSymbol):
1326    """
1327    Determines supersymmetry information from superspace group name; currently only for (3+1) superlattices
1328
1329    :param SGData: space group data structure as defined in SpcGroup above (see :ref:`SGData<SGData_table>`).
1330    :param SSymbol: superspace group symbol extension (string) defining modulation direction & generator info.
1331    :returns: (SSGError,SSGData)
1332   
1333       * SGError = 0 for no errors; >0 for errors (see SGErrors below for details)
1334       * SSGData - is a dict (see :ref:`Superspace Group object<SSGData_table>`) with entries:
1335       
1336             * 'SSpGrp': full superspace group symbol, accidental spaces removed; for display only
1337             * 'SSGCen': 4D cell centering vectors [0,0,0,0] at least
1338             * 'SSGOps': 4D symmetry operations as [M,T] so that M*x+T = x'
1339
1340    """
1341           
1342    def fixMonoOrtho():
1343        mod = ''.join(modsym).replace('1/2','0').replace('1','0')
1344        if SGData['SGPtGrp'] in ['2','m']:  #OK
1345            if mod in ['a00','0b0','00g']:
1346                result = [i*-1 for i in SGData['SSGKl']]
1347            else:
1348                result = SGData['SSGKl'][:]
1349            if '/' in mod:
1350                return [i*-1 for i in result]
1351            else:
1352                return result
1353        elif SGData['SGPtGrp'] == '2/m':    #OK
1354            if mod in ['a00','0b0','00g']:
1355                result =  SGData['SSGKl'][:]
1356            else:
1357                result = [i*-1 for i in SGData['SSGKl']]
1358            if '/' in mod:
1359                return [i*-1 for i in result]
1360            else:
1361                return result
1362        else:   #orthorhombic
1363            return [-SSGKl[i] if mod[i] in ['a','b','g'] else SSGKl[i] for i in range(3)]
1364       
1365    def extendSSGOps(SSGOps):
1366        for OpA in SSGOps:
1367            OpAtxt = SSMT2text(OpA)
1368            if 't' not in OpAtxt:
1369                continue
1370            for OpB in SSGOps:
1371                OpBtxt = SSMT2text(OpB)
1372                if 't' not in OpBtxt:
1373                    continue
1374                OpC = list(SGProd(OpB,OpA))
1375                OpC[1] %= 1.
1376                OpCtxt = SSMT2text(OpC)
1377#                print OpAtxt.replace(' ','')+' * '+OpBtxt.replace(' ','')+' = '+OpCtxt.replace(' ','')
1378                for k,OpD in enumerate(SSGOps):
1379                    OpDtxt = SSMT2text(OpD)
1380                    OpDtxt2 = ''
1381                    if SGData['SGGray']:                       
1382                        OpDtxt2 = SSMT2text([OpD[0],OpD[1]+np.array([0.,0.,0.,.5])])
1383#                    print '    ('+OpCtxt.replace(' ','')+' = ? '+OpDtxt.replace(' ','')+')'
1384                    if OpCtxt == OpDtxt:
1385                        continue
1386                    elif OpCtxt == OpDtxt2:
1387                        continue
1388                    elif OpCtxt.split(',')[:3] == OpDtxt.split(',')[:3]:
1389                        if 't' not in OpDtxt:
1390                            SSGOps[k] = OpC
1391#                            print k,'   new:',OpCtxt.replace(' ','')
1392                            break
1393                        else:
1394                            OpCtxt = OpCtxt.replace(' ','')
1395                            OpDtxt = OpDtxt.replace(' ','')
1396                            Txt = OpCtxt+' conflicts with '+OpDtxt
1397#                            print (Txt)
1398                            return False,Txt
1399        return True,SSGOps
1400       
1401    def findMod(modSym):
1402        for a in ['a','b','g']:
1403            if a in modSym:
1404                return a
1405               
1406    def genSSGOps():
1407        SSGOps = SSGData['SSGOps'][:]
1408        iFrac = {}
1409        for i,frac in enumerate(SSGData['modSymb']):
1410            if frac in ['1/2','1/3','1/4','1/6','1']:
1411                iFrac[i] = frac+'.'
1412#        print SGData['SpGrp']+SSymbol
1413#        print 'SSGKl',SSGKl,'genQ',genQ,'iFrac',iFrac,'modSymb',SSGData['modSymb']
1414# set identity & 1,-1; triclinic
1415        SSGOps[0][0][3,3] = 1.
1416## expand if centrosymmetric
1417#        if SGData['SGInv']:
1418#            SSGOps += [[-1*M,V] for M,V in SSGOps[:]]
1419# monoclinic - all done & all checked
1420        if SGData['SGPtGrp'] in ['2','m']:  #OK
1421            SSGOps[1][0][3,3] = SSGKl[0]
1422            SSGOps[1][1][3] = genQ[0]
1423            for i in iFrac:
1424                SSGOps[1][0][3,i] = -SSGKl[0]
1425        elif SGData['SGPtGrp'] == '2/m':    #OK
1426            SSGOps[1][0][3,3] = SSGKl[1]
1427            if gensym:
1428                SSGOps[1][1][3] = 0.5
1429            for i in iFrac:
1430                SSGOps[1][0][3,i] = SSGKl[0]
1431           
1432# orthorhombic - all OK not fully checked
1433        elif SGData['SGPtGrp'] in ['222','mm2','m2m','2mm']:    #OK
1434            if SGData['SGPtGrp'] == '222':
1435                OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[1,2],2:[1,3]},'b':{2:[3,2],0:[1,2]}} #OK
1436            elif SGData['SGPtGrp'] == 'mm2':
1437                OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]}} #OK
1438            elif SGData['SGPtGrp'] == 'm2m':
1439                OrOps = {'b':{0:[1,2],2:[3,2]},'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]}} #OK
1440            elif SGData['SGPtGrp'] == '2mm':
1441                OrOps = {'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]},'g':{0:[1,3],1:[2,3]}} #OK
1442            a = findMod(SSGData['modSymb'])
1443            OrFrac = OrOps[a]
1444            for j in iFrac:
1445                for i in OrFrac[j]:
1446                    SSGOps[i][0][3,j] = -2.*eval(iFrac[j])*SSGKl[i-1]
1447            for i in [0,1,2]:
1448                SSGOps[i+1][0][3,3] = SSGKl[i]
1449                SSGOps[i+1][1][3] = genQ[i]
1450                E,SSGOps = extendSSGOps(SSGOps)
1451                if not E:
1452                    return E,SSGOps
1453        elif SGData['SGPtGrp'] == 'mmm':    #OK
1454            OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]}} 
1455            a = findMod(SSGData['modSymb'])
1456            if a == 'g':
1457                SSkl = [1,1,1]
1458            elif a == 'a':
1459                SSkl = [-1,1,-1]
1460            else:
1461                SSkl = [1,-1,-1]
1462            OrFrac = OrOps[a]
1463            for j in iFrac:
1464                for i in OrFrac[j]:
1465                    SSGOps[i][0][3,j] = -2.*eval(iFrac[j])*SSkl[i-1]
1466            for i in [0,1,2]:
1467                SSGOps[i+1][0][3,3] = SSkl[i]
1468                SSGOps[i+1][1][3] = genQ[i]
1469                E,SSGOps = extendSSGOps(SSGOps)
1470                if not E:
1471                    return E,SSGOps               
1472# tetragonal - all done & checked
1473        elif SGData['SGPtGrp'] == '4':  #OK
1474            SSGOps[1][0][3,3] = SSGKl[0]
1475            SSGOps[1][1][3] = genQ[0]
1476            if '1/2' in SSGData['modSymb']:
1477                SSGOps[1][0][3,1] = -1
1478        elif SGData['SGPtGrp'] == '-4': #OK
1479            SSGOps[1][0][3,3] = SSGKl[0]
1480            if '1/2' in SSGData['modSymb']:
1481                SSGOps[1][0][3,1] = 1
1482        elif SGData['SGPtGrp'] in ['4/m',]: #OK
1483            if '1/2' in SSGData['modSymb']:
1484                SSGOps[1][0][3,1] = -SSGKl[0]
1485            for i,j in enumerate([1,3]):
1486                SSGOps[j][0][3,3] = 1
1487                if genQ[i]:
1488                    SSGOps[j][1][3] = genQ[i]
1489                E,SSGOps = extendSSGOps(SSGOps)
1490                if not E:
1491                    return E,SSGOps
1492        elif SGData['SGPtGrp'] in ['422','4mm','-42m','-4m2',]: #OK
1493            iGens = [1,4,5]
1494            if SGData['SGPtGrp'] in ['4mm','-4m2',]:
1495                iGens = [1,6,7]
1496            for i,j in enumerate(iGens):
1497                if '1/2' in SSGData['modSymb'] and i < 2:
1498                    SSGOps[j][0][3,1] = SSGKl[i]
1499                SSGOps[j][0][3,3] = SSGKl[i]
1500                if genQ[i]:
1501                    if 's' in gensym and j == 6:
1502                        SSGOps[j][1][3] = -genQ[i]
1503                    else:
1504                        SSGOps[j][1][3] = genQ[i]
1505                E,SSGOps = extendSSGOps(SSGOps)
1506                if not E:
1507                    return E,SSGOps
1508        elif SGData['SGPtGrp'] in ['4/mmm',]:#OK
1509            if '1/2' in SSGData['modSymb']:
1510                SSGOps[1][0][3,1] = -SSGKl[0]
1511                SSGOps[6][0][3,1] = SSGKl[1]
1512                if modsym:
1513                   SSGOps[1][1][3]  = -genQ[3]
1514            for i,j in enumerate([1,2,6,7]):
1515                SSGOps[j][0][3,3] = 1
1516                SSGOps[j][1][3] = genQ[i]
1517                E,Result = extendSSGOps(SSGOps)
1518                if not E:
1519                    return E,Result
1520                else:
1521                    SSGOps = Result
1522               
1523# trigonal - all done & checked
1524        elif SGData['SGPtGrp'] == '3':  #OK
1525            SSGOps[1][0][3,3] = SSGKl[0]
1526            if '1/3' in SSGData['modSymb']:
1527                SSGOps[1][0][3,1] = -1
1528            SSGOps[1][1][3] = genQ[0]
1529        elif SGData['SGPtGrp'] == '-3': #OK
1530            SSGOps[1][0][3,3] = -SSGKl[0]
1531            if '1/3' in SSGData['modSymb']:
1532                SSGOps[1][0][3,1] = -1
1533            SSGOps[1][1][3] = genQ[0]
1534        elif SGData['SGPtGrp'] in ['312','3m','-3m','-3m1','3m1']:   #OK
1535            if '1/3' in SSGData['modSymb']:
1536                SSGOps[1][0][3,1] = -1
1537            for i,j in enumerate([1,5]):
1538                if SGData['SGPtGrp'] in ['3m','-3m']:
1539                    SSGOps[j][0][3,3] = 1
1540                else:                   
1541                    SSGOps[j][0][3,3] = SSGKl[i+1]
1542                if genQ[i]:
1543                    SSGOps[j][1][3] = genQ[i]
1544        elif SGData['SGPtGrp'] in ['321','32']:   #OK
1545            for i,j in enumerate([1,4]):
1546                SSGOps[j][0][3,3] = SSGKl[i]
1547                if genQ[i]:
1548                    SSGOps[j][1][3] = genQ[i]
1549        elif SGData['SGPtGrp'] in ['31m','-31m']:   #OK
1550            ids = [1,3]
1551            if SGData['SGPtGrp'] == '-31m':
1552                ids = [1,3]
1553            if '1/3' in SSGData['modSymb']:
1554                SSGOps[ids[0]][0][3,1] = -SSGKl[0]
1555            for i,j in enumerate(ids):
1556                SSGOps[j][0][3,3] = 1
1557                if genQ[i+1]:
1558                    SSGOps[j][1][3] = genQ[i+1]
1559                     
1560# hexagonal all done & checked
1561        elif SGData['SGPtGrp'] == '6':  #OK
1562            SSGOps[1][0][3,3] = SSGKl[0]
1563            SSGOps[1][1][3] = genQ[0]
1564        elif SGData['SGPtGrp'] == '-6': #OK
1565            SSGOps[1][0][3,3] = SSGKl[0]
1566        elif SGData['SGPtGrp'] in ['6/m',]: #OK
1567            SSGOps[1][0][3,3] = -SSGKl[1]
1568            SSGOps[1][1][3] = genQ[0]
1569            SSGOps[2][1][3] = genQ[1]
1570        elif SGData['SGPtGrp'] in ['622',]: #OK
1571            for i,j in enumerate([1,8,9]):
1572                SSGOps[j][0][3,3] = SSGKl[i]
1573                if genQ[i]:
1574                    SSGOps[j][1][3] = genQ[i]
1575                E,SSGOps = extendSSGOps(SSGOps)
1576           
1577        elif SGData['SGPtGrp'] in ['6mm','-62m','-6m2',]: #OK
1578            for i,j in enumerate([1,6,7]):
1579                SSGOps[j][0][3,3] = SSGKl[i]
1580                if genQ[i]:
1581                    SSGOps[j][1][3] = genQ[i]
1582                E,SSGOps = extendSSGOps(SSGOps)
1583        elif SGData['SGPtGrp'] in ['6/mmm',]: # OK
1584            for i,j in enumerate([1,2,10,11]):
1585                SSGOps[j][0][3,3] = 1
1586                if genQ[i]:
1587                    SSGOps[j][1][3] = genQ[i]
1588                E,SSGOps = extendSSGOps(SSGOps)
1589        elif SGData['SGPtGrp'] in ['1','-1']: #triclinic - done
1590            return True,SSGOps
1591        E,SSGOps = extendSSGOps(SSGOps)
1592        return E,SSGOps
1593       
1594    def specialGen(gensym,modsym):
1595        sym = ''.join(gensym)
1596        if SGData['SGPtGrp'] in ['2/m',] and 'n' in SGData['SpGrp']:
1597            if 's' in sym:
1598                gensym = 'ss'
1599        if SGData['SGPtGrp'] in ['-62m',] and sym == '00s':
1600            gensym = '0ss'
1601        elif SGData['SGPtGrp'] in ['222',]:
1602            if sym == '00s':
1603                gensym = '0ss'
1604            elif sym == '0s0':
1605                gensym = 'ss0'
1606            elif sym == 's00':
1607                gensym = 's0s'
1608        elif SGData['SGPtGrp'] in ['mmm',]:
1609            if 'g' in modsym:
1610                if sym == 's00':
1611                    gensym = 's0s'
1612                elif sym == '0s0':
1613                    gensym = '0ss'
1614            elif 'a' in modsym:
1615                if sym == '0s0':
1616                    gensym = 'ss0'
1617                elif sym == '00s':
1618                    gensym = 's0s'
1619            elif 'b' in modsym:
1620                if sym == '00s':
1621                    gensym = '0ss'
1622                elif sym == 's00':
1623                    gensym = 'ss0'
1624        return gensym
1625                           
1626    Fracs = {'1/2':0.5,'1/3':1./3,'1':1.0,'0':0.,'s':.5,'t':1./3,'q':.25,'h':1./6,'a':0.,'b':0.,'g':0.}
1627    if SGData['SGLaue'] in ['m3','m3m']:
1628        return '(3+1) superlattices not defined for cubic space groups',None
1629    elif SGData['SGLaue'] in ['3R','3mR']:
1630        return '(3+1) superlattices not defined for rhombohedral settings - use hexagonal setting',None
1631    if SGData['SGGray'] and SSymbol[-1] == 's':
1632        SSymbol = SSymbol[:-1]
1633    try:
1634        modsym,gensym = splitSSsym(SSymbol)
1635    except ValueError:
1636        return 'Error in superspace symbol '+SSymbol,None
1637    modQ = [Fracs[mod] for mod in modsym]
1638    SSGKl = SGData['SSGKl'][:]
1639    if SGData['SGLaue'] in ['2/m','mmm']:
1640        SSGKl = fixMonoOrtho()
1641    Ngen = len(gensym)
1642    if len(gensym) and Ngen != len(SSGKl):
1643        return 'Wrong number of items in generator symbol '+''.join(gensym),None
1644    gensym = specialGen(gensym[:Ngen],modsym)
1645    genQ = [Fracs[mod] for mod in gensym[:Ngen]]
1646    if not genQ:
1647        genQ = [0,0,0,0]
1648    SSgSpc = SGData['SpGrp']+SSymbol
1649    if SGData['SGGray']:
1650        SSgSpc = SSgSpc.replace('('," 1'(")
1651        SSgSpc += 's'
1652    SSGData = {'SSpGrp':SSgSpc,'modQ':modQ,'modSymb':modsym,'SSGKl':SSGKl}
1653    SSCen = np.zeros((len(SGData['SGCen']),4))
1654    for icen,cen in enumerate(SGData['SGCen']):
1655        SSCen[icen,0:3] = cen
1656    SSCen[0] = np.zeros(4)
1657    SSGData['SSGCen'] = SSCen
1658    SSGData['SSGOps'] = []
1659    for iop,op in enumerate(SGData['SGOps']):
1660        T = np.zeros(4)
1661        ssop = np.zeros((4,4))
1662        ssop[:3,:3] = op[0]
1663        T[:3] = op[1]
1664        SSGData['SSGOps'].append([ssop,T])
1665    E,Result = genSSGOps()
1666    if E:
1667        SSGData['SSGOps'] = Result
1668        if DEBUG:
1669            print ('Super spacegroup operators for '+SSGData['SSpGrp'])
1670            for Op in Result:
1671                print (SSMT2text(Op).replace(' ',''))
1672            if SGData['SGInv']:                                 
1673                for Op in Result:
1674                    Op = [-Op[0],-Op[1]%1.]
1675                    print (SSMT2text(Op).replace(' ',''))                                 
1676        return None,SSGData
1677    else:
1678        return Result+'\nOperator conflict - incorrect superspace symbol',None
1679   
1680def SSChoice(SGData):
1681    '''
1682    Gets the unique set of possible super space groups for a given space group
1683    '''
1684    laueSS = {'-1':['(abg)',],
1685            '2/m':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1686            'mmm':['(00g)','(1/20g)','(01/2g)','(1/21/2g)','(10g)','(01g)',
1687                   '(a00)','(a1/20)','(a01/2)','(a1/21/2)','(a10)','(a01)',
1688                   '(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)','(1b0)','(0b1)',],
1689            '4/m':['(00g)','(1/21/2g)'],
1690            '4/mmm':['(00g)','(1/21/2g)'],
1691            '3':['(00g)','(1/31/3g)'],
1692            '3m1':['(00g)'],
1693            '31m':['(00g)','(1/31/3g)'],
1694            '6/m':['(00g)',],
1695            '6/mmm':['(00g)',]}
1696           
1697    laueTS = {'-1':['',],
1698            '2/m':['','s','0s','ss','s0'],
1699            'mmm':['','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1700            '4/m':['','q','s','s0',],
1701            '4/mmm':['','q00','s00','s0s','ss0','0ss','qq0','qqs','0q0','s0s0','00ss','s00s','ss00','0ss0','0s0s'],
1702            '3':['','t'],
1703            '3m1':['','t0','0s','t00','0s0'],
1704            '31m':['','t00','0ss'],
1705            '6/m':['','h','t','s','s0'],
1706            '6/mmm':['','h00','t00','s00','ss0','0ss','s0s','s0s0','00ss','s00s','ss00','0ss0','0s0s']}
1707    laue = SGData['SGLaue']
1708    SSChoice = []
1709    for ax in laueSS[laue]:
1710        for sx in laueTS[laue]:
1711            SSChoice.append(ax+sx)               
1712    ssChoice = []
1713    ssHash = []
1714    for item in SSChoice:
1715        E,SSG = SSpcGroup(SGData,item)
1716        if SSG:
1717            sshash = hash(str(SSGPrint(SGData,SSG)[1]))
1718            if sshash not in ssHash:
1719                ssHash.append(sshash)
1720                ssChoice.append(item)
1721    return ssChoice
1722
1723def fixGray(SGData,SSymbol):
1724    modsym,gensym = SSymbol.replace(' ','').split(')')
1725    modsym += ')'
1726    sgPtGp = SGData['SGPtGrp']
1727    if gensym:
1728        if sgPtGp in ['1','2','m','3','4','6'] and len(gensym) == 1:
1729            gensym += 's'
1730        elif sgPtGp in ['2/m','4/m','6/m'] and len(gensym) == 2:
1731            gensym += 's'
1732        elif sgPtGp in ['4/mmm','6/mmm'] and len(gensym) == 4:
1733            gensym += 's'
1734        elif len(gensym) == 3:
1735            gensym += 's'
1736    else:
1737        if sgPtGp in ['1','2','m','3','4','6']:
1738            gensym += '0s'
1739        elif sgPtGp in ['2/m','4/m','6/m']:
1740            gensym += '00s'
1741        elif sgPtGp in ['4/mmm','6/mmm']:
1742            gensym += '0000s'
1743        else:
1744            gensym += '000s'
1745    return modsym+gensym
1746           
1747def splitSSsym(SSymbol):
1748    '''
1749    Splits supersymmetry symbol into two lists of strings
1750    '''
1751    modsym,gensym = SSymbol.replace(' ','').split(')')
1752    modsym = modsym.replace(',','')
1753    if "1'" in modsym:
1754        gensym = gensym[:-1]
1755    modsym = modsym.replace("1'",'')
1756    if gensym in ['0','00','000','0000']:       #get rid of extraneous symbols
1757        gensym = ''
1758    nfrac = modsym.count('/')
1759    modsym = modsym.lstrip('(')
1760    if nfrac == 0:
1761        modsym = list(modsym)
1762    elif nfrac == 1:
1763        pos = modsym.find('/')
1764        if pos == 1:
1765            modsym = [modsym[:3],modsym[3],modsym[4]]
1766        elif pos == 2:
1767            modsym = [modsym[0],modsym[1:4],modsym[4]]
1768        else:
1769            modsym = [modsym[0],modsym[1],modsym[2:]]
1770    else:
1771        lpos = modsym.find('/')
1772        rpos = modsym.rfind('/')
1773        if lpos == 1 and rpos == 4:
1774            modsym = [modsym[:3],modsym[3:6],modsym[6]]
1775        elif lpos == 1 and rpos == 5:
1776            modsym = [modsym[:3],modsym[3],modsym[4:]]
1777        else:
1778            modsym = [modsym[0],modsym[1:4],modsym[4:]]
1779    gensym = list(gensym)
1780    return modsym,gensym
1781       
1782def SSGPrint(SGData,SSGData):
1783    '''
1784    Print the output of SSpcGroup in a nicely formatted way. Used in SSpaceGroup
1785
1786    :param SGData: space group data structure as defined in SpcGroup above.
1787    :param SSGData: from :func:`SSpcGroup`
1788    :returns:
1789        SSGText - list of strings with the superspace group details
1790        SGTable - list of strings for each of the operations
1791    '''
1792    Mult = len(SSGData['SSGCen'])*len(SSGData['SSGOps'])*(int(SGData['SGInv'])+1)
1793    SSGText = []
1794    SSGText.append(' Superspace Group: '+SSGData['SSpGrp'])
1795    CentStr = 'centrosymmetric'
1796    if not SGData['SGInv']:
1797        CentStr = 'non'+CentStr
1798    if SGData['SGLatt'] in 'ABCIFR':
1799        SSGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
1800    else:
1801        SSGText.append(' The superlattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())
1802    SSGText.append(' The Laue symmetry is '+SGData['SGLaue'])
1803    SGptGp = SGData['SGPtGrp']
1804    if SGData['SGGray']:
1805        SGptGp += "1'"
1806    SSGText.append(' The superlattice point group is '+SGptGp+', '+''.join([str(i) for i in SSGData['SSGKl']]))
1807    SSGText.append(' The number of superspace group generators is '+str(len(SGData['SSGKl'])))
1808    SSGText.append(' Multiplicity of a general site is '+str(Mult))
1809    if SGData['SGUniq'] in ['a','b','c']:
1810        SSGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
1811    if SGData['SGInv']:
1812        SSGText.append(' The inversion center is located at 0,0,0')
1813    if SGData['SGPolax']:
1814        SSGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
1815    SSGText.append(' ')
1816    if len(SSGData['SSGCen']) > 1:
1817        SSGText.append(' The equivalent positions are:')
1818        SSGText.append(' ('+SSLatt2text(SSGData['SSGCen'])+')+\n')
1819    else:
1820        SSGText.append(' The equivalent positions are:\n')
1821    SSGTable = []
1822    for i,Opr in enumerate(SSGData['SSGOps']):
1823        SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr)))
1824    if SGData['SGGray']:
1825        SSGTable.append("     for 1'")
1826        for i,Opr in enumerate(SSGData['SSGOps']):
1827            Opr2 = [Opr[0],Opr[1]+np.array([0,0,0,.5])]
1828            SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr2)))
1829    return SSGText,SSGTable
1830   
1831def SSGModCheck(Vec,modSymb,newMod=True):
1832    ''' Checks modulation vector compatibility with supersymmetry space group symbol.
1833    if newMod: Superspace group symbol takes precidence & the vector will be modified accordingly
1834    '''
1835    Fracs = {'1/2':0.5,'1/3':1./3,'1':1.0,'0':0.,'a':0.,'b':0.,'g':0.}
1836    modQ = [Fracs[mod] for mod in modSymb]
1837    if newMod:
1838        newVec = Vec
1839        if not np.any(Vec):
1840            newVec = [0.1 if (vec == 0.0 and mod in ['a','b','g']) else vec for [vec,mod] in zip(Vec,modSymb)]
1841        return [Q if mod not in ['a','b','g'] and vec != Q else vec for [vec,mod,Q] in zip(newVec,modSymb,modQ)],  \
1842            [True if mod in ['a','b','g'] else False for mod in modSymb]
1843    else:
1844        return Vec,[True if mod in ['a','b','g'] else False for mod in modSymb]
1845
1846def SSMT2text(Opr):
1847    "From superspace group matrix/translation operator returns text version"
1848    XYZS = ('x','y','z','t')    #Stokes, Campbell & van Smaalen notation
1849    TRA = ('   ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
1850    Fld = ''
1851    M,T = Opr
1852    for j in range(4):
1853        IJ = ''
1854        for k in range(4):
1855            txt = str(int(round(M[j][k])))
1856            txt = txt.replace('1',XYZS[k]).replace('0','')
1857            if '2' in txt:
1858                txt += XYZS[k]
1859            if IJ and M[j][k] > 0:
1860                IJ += '+'+txt
1861            else:
1862                IJ += txt
1863        IK = int(round(T[j]*12))%12
1864        if IK:
1865            if not IJ:
1866                break
1867            if IJ[0] == '-':
1868                Fld += (TRA[IK]+IJ).rjust(8)
1869            else:
1870                Fld += (TRA[IK]+'+'+IJ).rjust(8)
1871        else:
1872            Fld += IJ.rjust(8)
1873        if j != 3: Fld += ', '
1874    return Fld
1875   
1876def SSLatt2text(SSGCen):
1877    "Lattice centering vectors to text"
1878    lattTxt = ''
1879    lattDir = {4:'1/3',6:'1/2',8:'2/3',0:'0'}
1880    for vec in SSGCen:
1881        lattTxt += ' '
1882        for item in vec:
1883            lattTxt += '%s,'%(lattDir[int(item*12)])
1884        lattTxt = lattTxt.rstrip(',')
1885        lattTxt += ';'
1886    lattTxt = lattTxt.rstrip(';').lstrip(' ')
1887    return lattTxt
1888       
1889def SSpaceGroup(SGSymbol,SSymbol):
1890    '''
1891    Print the output of SSpcGroup in a nicely formatted way.
1892
1893    :param SGSymbol: space group symbol with spaces between axial fields.
1894    :param SSymbol: superspace group symbol extension (string).
1895    :returns: nothing
1896    '''
1897
1898    E,A = SpcGroup(SGSymbol)
1899    if E > 0:
1900        print (SGErrors(E))
1901        return
1902    E,B = SSpcGroup(A,SSymbol)   
1903    if E > 0:
1904        print (E)
1905        return
1906    for l in SSGPrint(B):
1907        print (l)
1908       
1909def SGProd(OpA,OpB):
1910    '''
1911    Form space group operator product. OpA & OpB are [M,V] pairs;
1912        both must be of same dimension (3 or 4). Returns [M,V] pair
1913    '''
1914    A,U = OpA
1915    B,V = OpB
1916    M = np.inner(B,A.T)
1917    W = np.inner(B,U)+V
1918    return M,W
1919       
1920def GetLittleGrpOps(SGData,vec):
1921    ''' Find rotation part of operators that leave vec unchanged
1922   
1923    :param SGData: space group data structure as defined in SpcGroup above.
1924    :param vec: a numpy array of fractional vector coordinates
1925    :returns: Little - list of operators [M,T] that form the little gropu
1926    '''
1927    Little = []
1928    Ops = SGData['SGOps'][:]
1929    if SGData['SGInv']:
1930        Ops += [[-M,-T] for [M,T] in Ops]
1931    for [M,T] in Ops:
1932        tvec = np.inner(M,vec)%1.
1933        if np.allclose(tvec,vec%1.):
1934            Little.append([M,T])
1935    return Little
1936       
1937def MoveToUnitCell(xyz):
1938    '''
1939    Translates a set of coordinates so that all values are >=0 and < 1
1940
1941    :param xyz: a list or numpy array of fractional coordinates
1942    :returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
1943    '''
1944    XYZ = (np.array(xyz)+10.)%1.
1945    cell = np.asarray(np.rint(XYZ-xyz),dtype=np.int32)
1946    return XYZ,cell
1947       
1948def Opposite(XYZ,toler=0.0002):
1949    '''
1950    Gives opposite corner, edge or face of unit cell for position within tolerance.
1951        Result may be just outside the cell within tolerance
1952
1953    :param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
1954    :param toler: unit cell fraction tolerance making opposite
1955    :returns:
1956        XYZ: dict of opposite positions; key=unit cell & always contains XYZ
1957    '''
1958    perm3 = [[1,1,1],[0,1,1],[1,0,1],[1,1,0],[1,0,0],[0,1,0],[0,0,1],[0,0,0]]
1959    TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
1960    perm = TB*perm3
1961    cperm = ['%d,%d,%d'%(i,j,k) for i,j,k in perm]
1962    D = dict(zip(cperm,perm))
1963    new = {}
1964    for key in D:
1965        new[key] = np.array(D[key])+np.array(XYZ)
1966    return new
1967       
1968def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
1969    '''
1970    Generates the equivalent positions for a specified coordinate and space group
1971
1972    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
1973    :param SGData: from :func:`SpcGroup`
1974    :param All: True return all equivalent positions including duplicates;
1975      False return only unique positions
1976    :param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
1977    :param Move: True move generated atom positions to be inside cell
1978      False do not move atoms       
1979    :return: [[XYZEquiv],Idup,[UijEquiv],spnflp]
1980
1981      *  [XYZEquiv] is list of equivalent positions (XYZ is first entry)
1982      *  Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
1983      * is centering operator number (1-4) and - is for inversion
1984        Cell = unit cell translations needed to put new positions inside cell
1985        [UijEquiv] - equivalent Uij; absent if no Uij given
1986      * +1/-1 for spin inversion of operator - empty if not magnetic
1987       
1988    '''
1989    XYZEquiv = []
1990    UijEquiv = []
1991    Idup = []
1992    Cell = []
1993    inv = int(SGData['SGInv']+1)
1994    icen = SGData['SGCen']
1995    if SGData['SGFixed']:
1996        inv = 1
1997    SpnFlp = SGData.get('SpnFlp',[])
1998    spnflp = []
1999    X = np.array(XYZ)
2000    mj = 0
2001    cell0 = np.zeros(3,dtype=np.int32)
2002    if Move:
2003        X,cell0 = MoveToUnitCell(X)
2004    for ic,cen in enumerate(icen):
2005        C = np.array(cen)
2006        for invers in range(inv):
2007            for io,[M,T] in enumerate(SGData['SGOps']):
2008                idup = ((io+1)+100*ic)*(1-2*invers)
2009                XT = np.inner(M,X)+T
2010                if len(Uij):
2011                    U = Uij2U(Uij)
2012                    U = np.inner(M,np.inner(U,M).T)
2013                    newUij = U2Uij(U)
2014                if invers:
2015                    XT = -XT
2016                XT += C
2017                cell = np.zeros(3,dtype=np.int32)+cell0
2018                cellj = np.zeros(3,dtype=np.int32)
2019                if Move:
2020                    newX,cellj = MoveToUnitCell(XT)
2021                else:
2022                    newX = XT
2023                cell += cellj
2024                if All:
2025                    if np.allclose(newX,X,atol=0.0002):
2026                        idup = False
2027                else:
2028                    if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
2029                        idup = False
2030                if All or idup:
2031                    XYZEquiv.append(newX)
2032                    Idup.append(idup)
2033                    Cell.append(cell)
2034                    if len(Uij):
2035                        UijEquiv.append(newUij)
2036                    if len(SpnFlp):
2037                        spnflp.append(SpnFlp[mj])
2038                    else:
2039                        spnflp.append(1)
2040                mj += 1
2041    if len(Uij):
2042        return zip(XYZEquiv,UijEquiv,Idup,Cell,spnflp)
2043    else:
2044        return zip(XYZEquiv,Idup,Cell,spnflp)
2045       
2046def GenHKL(HKL,SGData):
2047    ''' Generates all equivlent reflections including Friedel pairs
2048    :param HKL:  [h,k,l] must be integral values
2049    :param SGData: space group data obtained from SpcGroup
2050    :returns: array Uniq: equivalent reflections
2051    '''
2052   
2053    Ops = SGData['SGOps']
2054    OpM = np.array([op[0] for op in Ops])
2055    Uniq = np.inner(OpM,HKL)
2056    Uniq = list(Uniq)+list(-1*Uniq)
2057    return np.array(Uniq)
2058
2059def GenHKLf(HKL,SGData):
2060    '''
2061    Uses old GSAS Fortran routine genhkl.for
2062
2063    :param HKL:  [h,k,l] must be integral values for genhkl.for to work
2064    :param SGData: space group data obtained from SpcGroup
2065    :returns: iabsnt,mulp,Uniq,phi
2066
2067     *   iabsnt = True if reflection is forbidden by symmetry
2068     *   mulp = reflection multiplicity including Friedel pairs
2069     *   Uniq = numpy array of equivalent hkl in descending order of h,k,l
2070     *   phi = phase offset for each equivalent h,k,l
2071
2072    '''
2073    hklf = list(HKL)+[0,]       #could be numpy array!
2074    Ops = SGData['SGOps']
2075    OpM = np.array([op[0] for op in Ops],order='F')
2076    OpT = np.array([op[1] for op in Ops])
2077    Cen = np.array([cen for cen in SGData['SGCen']],order='F')
2078   
2079    import pyspg
2080    Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
2081    h,k,l,f = Uniq
2082    Uniq=np.array(list(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq])))
2083    phi = f[:Nuniq]
2084    return iabsnt,mulp,Uniq,phi
2085   
2086def checkSSLaue(HKL,SGData,SSGData):
2087    #Laue check here - Toss HKL if outside unique Laue part
2088    h,k,l,m = HKL
2089    if SGData['SGLaue'] == '2/m':
2090        if SGData['SGUniq'] == 'a':
2091            if 'a' in SSGData['modSymb'] and h == 0 and m < 0:
2092                return False
2093            elif 'b' in SSGData['modSymb'] and k == 0 and l ==0 and m < 0:
2094                return False
2095            else:
2096                return True
2097        elif SGData['SGUniq'] == 'b':
2098            if 'b' in SSGData['modSymb'] and k == 0 and m < 0:
2099                return False
2100            elif 'a' in SSGData['modSymb'] and h == 0 and l ==0 and m < 0:
2101                return False
2102            else:
2103                return True
2104        elif SGData['SGUniq'] == 'c':
2105            if 'g' in SSGData['modSymb'] and l == 0 and m < 0:
2106                return False
2107            elif 'a' in SSGData['modSymb'] and h == 0 and k ==0 and m < 0:
2108                return False
2109            else:
2110                return True
2111    elif SGData['SGLaue'] == 'mmm':
2112        if 'a' in SSGData['modSymb']:
2113            if h == 0 and m < 0:
2114                return False
2115            else:
2116                return True
2117        elif 'b' in SSGData['modSymb']:
2118            if k == 0 and m < 0:
2119                return False
2120            else:
2121                return True
2122        elif 'g' in SSGData['modSymb']:
2123            if l == 0 and m < 0:
2124                return False
2125            else:
2126                return True
2127    else:   #tetragonal, trigonal, hexagonal (& triclinic?)
2128        if l == 0 and m < 0:
2129            return False
2130        else:
2131            return True
2132       
2133def checkHKLextc(HKL,SGData):
2134    '''
2135    Checks if reflection extinct - does not check centering
2136
2137    :param HKL:  [h,k,l]
2138    :param SGData: space group data obtained from SpcGroup
2139    :returns: True if extinct; False if allowed
2140
2141    '''
2142    Ops = SGData['SGOps']
2143    OpM = np.array([op[0] for op in Ops])
2144    OpT = np.array([op[1] for op in Ops])
2145    HKLS = np.array([HKL,-HKL])     #Freidel's Law
2146    DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,3))
2147    PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2148    for dhkl,phkl in zip(DHKL,PHKL)[1:]:    #skip identity
2149        if dhkl.any():
2150            continue
2151        else:
2152            if phkl%1.:
2153                return True
2154    return False
2155
2156def checkMagextc(HKL,SGData):
2157    '''
2158    Checks if reflection magnetically extinct; does fullcheck (centering, too)
2159    uses algorthm from Gallego, et al., J. Appl. Cryst. 45, 1236-1247 (2012)
2160
2161    :param HKL:  [h,k,l]
2162    :param SGData: space group data obtained from SpcGroup; must have magnetic symmetry SpnFlp data
2163    :returns: True if magnetically extinct; False if allowed (to match GenHKLf)
2164
2165    '''
2166    Ops = SGData['SGOps']
2167    Ncen = len(SGData['SGCen'])
2168    OpM = np.array([op[0] for op in Ops])
2169    OpT = np.array([op[1] for op in Ops])
2170    if SGData['SGInv'] and not SGData['SGFixed']:
2171        OpM = np.vstack((OpM,-OpM))
2172        OpT = np.vstack((OpT,-OpT))%1.
2173    OpM = np.reshape(np.array(list(OpM)*Ncen),(-1,3,3))
2174    OpT = np.reshape(np.array([OpT+cen for cen in SGData['SGCen']]),(-1,3))
2175    Spn = SGData['SpnFlp'][:len(OpM)]
2176    Mag = np.array([nl.det(opm) for opm in OpM])*Spn
2177    DHKL = np.reshape(np.inner(HKL,OpM),(-1,3))
2178    PHKL = np.reshape(np.cos(twopi*np.inner(HKL,OpT))*Mag,(-1,))[:,nxs,nxs]*OpM     #compute T(R,theta) eq(7)
2179    Ftest = np.random.rand(3)       #random magnetic moment
2180    Psum = np.zeros(3)
2181    nsum = 0.
2182    nA = 0
2183    for dhkl,phkl in zip(DHKL,PHKL):
2184        if not np.allclose(dhkl,HKL):           #test for eq(5)
2185            continue
2186        else:
2187            nA += 1
2188            nsum += np.trace(phkl)          #eq(8)
2189            pterm = np.inner(Ftest,phkl)    #eq(9)
2190            Psum += pterm
2191    if nsum/nA > 1.:        #only need to look at nA=1 frok eq(8)
2192        return False
2193    if np.allclose(Psum,np.zeros(3)):
2194        return True
2195    else:
2196        if np.inner(HKL,Psum):
2197            return True
2198        return False
2199   
2200def checkSSextc(HKL,SSGData):
2201    Ops = SSGData['SSGOps']
2202    OpM = np.array([op[0] for op in Ops])
2203    OpT = np.array([op[1] for op in Ops])
2204    HKLS = np.array([HKL,-HKL])     #Freidel's Law
2205    DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,4))
2206    PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2207    for dhkl,phkl in zip(DHKL,PHKL)[1:]:    #skip identity
2208        if dhkl.any():
2209            continue
2210        else:
2211            if phkl%1.:
2212                return False
2213    return True
2214   
2215################################################################################
2216#### Site symmetry tables
2217################################################################################
2218     
2219OprName = {
2220    '-6643':       ['-1',1],'6479' :    ['2(z)',2],'-6479':     ['m(z)',3],
2221    '6481' :     ['m(y)',4],'-6481':    ['2(y)',5],'6641' :     ['m(x)',6],
2222    '-6641':     ['2(x)',7],'6591' :  ['m(+-0)',8],'-6591':   ['2(+-0)',9],
2223    '6531' :  ['m(110)',10],'-6531': ['2(110)',11],'6537' :    ['4(z)',12],
2224    '-6537':   ['-4(z)',13],'975'  : ['3(111)',14],'6456' :       ['3',15],
2225    '-489' :  ['3(+--)',16],'483'  : ['3(-+-)',17],'-969' :  ['3(--+)',18],
2226    '819'  :  ['m(+0-)',19],'-819' : ['2(+0-)',20],'2431' :  ['m(0+-)',21],
2227    '-2431':  ['2(0+-)',22],'-657' :  ['m(xz)',23],'657'  :   ['2(xz)',24],
2228    '1943' :   ['-4(x)',25],'-1943':   ['4(x)',26],'-2429':   ['m(yz)',27],
2229    '2429' :   ['2(yz)',28],'639'  :  ['-4(y)',29],'-639' :    ['4(y)',30],
2230    '-6484':   ['2(010)',4],'6484' :  ['m(010)',5],'-6668':   ['2(100)',6],
2231    '6668' :   ['m(100)',7],'-6454': ['2(120)',18],'6454' :  ['m(120)',19],
2232    '-6638':  ['2(210)',20],'6638' : ['m(210)',21],   #search in SytSym ends at m(210)
2233    '2223' : ['3(+++)2',39],
2234    '6538' :   ['6(z)1',40],'-2169':['3(--+)2',41],'2151' : ['3(+--)2',42],
2235    '2205' :['-3(-+-)2',43],'-2205':[' (-+-)2',44],'489'  :['-3(+--)1',45],
2236    '801'  :   ['4(y)1',46],'1945' :  ['4(x)3',47],'-6585': ['-4(z)3 ',48],
2237    '6585' :   ['4(z)3',49],'6584' :  ['3(z)2',50],'6666' :  ['6(z)5 ',51],
2238    '6643' :       ['1',52],'-801' : ['-4(y)1',53],'-1945': ['-4(x)3 ',54],
2239    '-6666':  ['-6(z)5',55],'-6538': ['-6(z)1',56],'-2223':['-3(+++)2',57],
2240    '-975' :['-3(+++)1',58],'-6456': ['-3(z)1',59],'-483' :['-3(-+-)1',60],
2241    '969'  :['-3(--+)1',61],'-6584': ['-3(z)2',62],'2169' :['-3(--+)2',63],
2242    '-2151':['-3(+--)2',64],   }                               
2243
2244KNsym = {
2245    '0'         :'    1   ','1'         :'   -1   ','64'        :'    2(x)','32'        :'    m(x)',
2246    '97'        :'  2/m(x)','16'        :'    2(y)','8'         :'    m(y)','25'        :'  2/m(y)',
2247    '2'         :'    2(z)','4'         :'    m(z)','7'         :'  2/m(z)','134217728' :'   2(yz)',
2248    '67108864'  :'   m(yz)','201326593' :' 2/m(yz)','2097152'   :'  2(0+-)','1048576'   :'  m(0+-)',
2249    '3145729'   :'2/m(0+-)','8388608'   :'   2(xz)','4194304'   :'   m(xz)','12582913'  :' 2/m(xz)',
2250    '524288'    :'  2(+0-)','262144'    :'  m(+0-)','796433'    :'2/m(+0-)','1024'      :'   2(xy)',
2251    '512'       :'   m(xy)','1537'      :' 2/m(xy)','256'       :'  2(+-0)','128'       :'  m(+-0)',
2252    '385'       :'2/m(+-0)','76'        :'  mm2(x)','52'        :'  mm2(y)','42'        :'  mm2(z)',
2253    '135266336' :' mm2(yz)','69206048'  :'mm2(0+-)','8650760'   :' mm2(xz)','4718600'   :'mm2(+0-)',
2254    '1156'      :' mm2(xy)','772'       :'mm2(+-0)','82'        :'  222   ','136314944' :'  222(x)',
2255    '8912912'   :'  222(y)','1282'      :'  222(z)','127'       :'  mmm   ','204472417' :'  mmm(x)',
2256    '13369369'  :'  mmm(y)','1927'      :'  mmm(z)','33554496'  :'  4(x)','16777280'  :' -4(x)',
2257    '50331745'  :'4/m(x)'  ,'169869394' :'422(x)','84934738'  :'-42m(x)','101711948' :'4mm(x)',
2258    '254804095' :'4/mmm(x)','536870928 ':'  4(y)','268435472' :' -4(y)','805306393' :'4/m(y)',
2259    '545783890' :'422(y)','272891986' :'-42m(y)','541327412' :'4mm(y)','818675839' :'4/mmm(y)',
2260    '2050'      :'  4(z)','4098'      :' -4(z)','6151'      :'4/m(z)','3410'      :'422(z)',
2261    '4818'      :'-42m(z)','2730'      :'4mm(z)','8191'      :'4/mmm(z)','8192'      :'  3(111)',
2262    '8193'      :' -3(111)','2629888'   :' 32(111)','1319040'   :' 3m(111)','3940737'   :'-3m(111)',
2263    '32768'     :'  3(+--)','32769'     :' -3(+--)','10519552'  :' 32(+--)','5276160'   :' 3m(+--)',
2264    '15762945'  :'-3m(+--)','65536'     :'  3(-+-)','65537'     :' -3(-+-)','134808576' :' 32(-+-)',
2265    '67437056'  :' 3m(-+-)','202180097' :'-3m(-+-)','131072'    :'  3(--+)','131073'    :' -3(--+)',
2266    '142737664' :' 32(--+)','71434368'  :' 3m(--+)','214040961' :'-3m(--+)','237650'    :'   23   ',
2267    '237695'    :'   m3   ','715894098' :'   432  ','358068946' :'  -43m  ','1073725439':'   m3m  ',
2268    '68157504'  :' mm2(d100)','4456464'   :' mm2(d010)','642'       :' mm2(d001)','153092172' :'-4m2(x)',
2269    '277348404' :'-4m2(y)','5418'      :'-4m2(z)','1075726335':'  6/mmm ','1074414420':'-6m2(100)',
2270    '1075070124':'-6m2(120)','1075069650':'   6mm  ','1074414890':'   622  ','1073758215':'   6/m  ',
2271    '1073758212':'   -6   ','1073758210':'    6   ','1073759865':'-3m(100)','1075724673':'-3m(120)',
2272    '1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
2273    '1073758209':'   -3   ','1073758208':'    3   ','1074135143':'mmm(100)','1075314719':'mmm(010)',
2274    '1073743751':'mmm(110)','1074004034':' mm2(z100)','1074790418':' mm2(z010)','1073742466':' mm2(z110)',
2275    '1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
2276    '1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
2277    '1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
2278    '1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
2279    '1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
2280    '1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
2281    '1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
2282    '1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':'   -1   ',
2283    }
2284
2285NXUPQsym = {
2286    '1'        :(28,29,28,28),'-1'       :( 1,29,28, 0),'2(x)'     :(12,18,12,25),'m(x)'     :(25,18,12,25),
2287    '2/m(x)'   :( 1,18, 0,-1),'2(y)'     :(13,17,13,24),'m(y)'     :(24,17,13,24),'2/m(y)'   :( 1,17, 0,-1),
2288    '2(z)'     :(14,16,14,23),'m(z)'     :(23,16,14,23),'2/m(z)'   :( 1,16, 0,-1),'2(yz)'    :(10,23,10,22),
2289    'm(yz)'    :(22,23,10,22),' 2/m(yz)' :( 1,23, 0,-1),'2(0+-)'   :(11,24,11,21),'m(0+-)'   :(21,24,11,21),
2290    '2/m(0+-)' :( 1,24, 0,-1),'2(xz)'    :( 8,21, 8,20),'m(xz)'    :(20,21, 8,20),'2/m(xz)'  :( 1,21, 0,-1),
2291    '2(+0-)'   :( 9,22, 9,19),'m(+0-)'   :(19,22, 9,19),'2/m(+0-)' :( 1,22, 0,-1),'2(xy)'    :( 6,19, 6,18),
2292    'm(xy)'    :(18,19, 6,18),' 2/m(xy)' :( 1,19, 0,-1),'2(+-0)'   :( 7,20, 7,17),'m(+-0)'   :(17,20, 7,17),
2293    '2/m(+-0)' :( 1,20, 17,-1),'mm2(x)'  :(12,10, 0,-1),'mm2(y)'   :(13,10, 0,-1),'mm2(z)'   :(14,10, 0,-1),
2294    'mm2(yz)'  :(10,13, 0,-1),'mm2(0+-)' :(11,13, 0,-1),'mm2(xz)'  :( 8,12, 0,-1),'mm2(+0-)' :( 9,12, 0,-1),
2295    'mm2(xy)'  :( 6,11, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222'      :( 1,10, 0,-1),'222(x)'   :( 1,13, 0,-1),
2296    '222(y)'   :( 1,12, 0,-1),'222(z)'   :( 1,11, 0,-1),'mmm'      :( 1,10, 0,-1),'mmm(x)'   :( 1,13, 0,-1),
2297    'mmm(y)'   :( 1,12, 0,-1),'mmm(z)'   :( 1,11, 0,-1),'4(x)'     :(12, 4,12, 0),'-4(x)'    :( 1, 4,12, 0),
2298    '4/m(x)'   :( 1, 4,12,-1),'422(x)'   :( 1, 4, 0,-1),'-42m(x)'  :( 1, 4, 0,-1),'4mm(x)'   :(12, 4, 0,-1),
2299    '4/mmm(x)' :( 1, 4, 0,-1),'4(y)'     :(13, 3,13, 0),'-4(y)'    :( 1, 3,13, 0),'4/m(y)'   :( 1, 3,13,-1),
2300    '422(y)'   :( 1, 3, 0,-1),'-42m(y)'  :( 1, 3, 0,-1),'4mm(y)'   :(13, 3, 0,-1),'4/mmm(y)' :(1, 3, 0,-1,),
2301    '4(z)'     :(14, 2,14, 0),'-4(z)'    :( 1, 2,14, 0),'4/m(z)'   :( 1, 2,14,-1),'422(z)'   :( 1, 2, 0,-1),
2302    '-42m(z)'  :( 1, 2, 0,-1),'4mm(z)'   :(14, 2, 0,-1),'4/mmm(z)' :( 1, 2, 0,-1),'3(111)'   :( 2, 5, 2, 0),
2303    '-3(111)'  :( 1, 5, 2, 0),'32(111)'  :( 1, 5, 0, 2),'3m(111)'  :( 2, 5, 0, 2),'-3m(111)' :( 1, 5, 0,-1),
2304    '3(+--)'   :( 5, 8, 5, 0),'-3(+--)'  :( 1, 8, 5, 0),'32(+--)'  :( 1, 8, 0, 5),'3m(+--)'  :( 5, 8, 0, 5),
2305    '-3m(+--)' :( 1, 8, 0,-1),'3(-+-)'   :( 4, 7, 4, 0),'-3(-+-)'  :( 1, 7, 4, 0),'32(-+-)'  :( 1, 7, 0, 4),
2306    '3m(-+-)'  :( 4, 7, 0, 4),'-3m(-+-)' :( 1, 7, 0,-1),'3(--+)'   :( 3, 6, 3, 0),'-3(--+)'  :( 1, 6, 3, 0),
2307    '32(--+)'  :( 1, 6, 0, 3),'3m(--+)'  :( 3, 6, 0, 3),'-3m(--+)' :( 1, 6, 0,-1),'23'       :( 1, 1, 0, 0),
2308    'm3'       :( 1, 1, 0, 0),'432'      :( 1, 1, 0, 0),'-43m'     :( 1, 1, 0, 0),'m3m'      :( 1, 1, 0, 0),
2309    'mm2(d100)':(12,13, 0,-1),'mm2(d010)':(13,12, 0,-1),'mm2(d001)':(14,11, 0,-1),'-4m2(x)'  :( 1, 4, 0,-1),
2310    '-4m2(y)'  :( 1, 3, 0,-1),'-4m2(z)'  :( 1, 2, 0,-1),'6/mmm'    :( 1, 9, 0,-1),'-6m2(100)':( 1, 9, 0,-1),
2311    '-6m2(120)':( 1, 9, 0,-1),'6mm'      :(14, 9, 0,-1),'622'      :( 1, 9, 0,-1),'6/m'      :( 1, 9,14,-1),
2312    '-6'       :( 1, 9,14, 0),'6'        :(14, 9,14, 0),'-3m(100)' :( 1, 9, 0,-1),'-3m(120)' :( 1, 9, 0,-1),
2313    '3m(100)'  :(14, 9, 0,14),'3m(120)'  :(14, 9, 0,14),'32(100)'  :( 1, 9, 0,14),'32(120)'  :( 1, 9, 0,14),
2314    '-3'       :( 1, 9,14, 0),'3'        :(14, 9,14, 0),'mmm(100)' :( 1,14, 0,-1),'mmm(010)' :( 1,15, 0,-1),
2315    'mmm(110)' :( 1,11, 0,-1),'mm2(z100)':(14,14, 0,-1),'mm2(z010)':(14,15, 0,-1),'mm2(z110)':(14,11, 0,-1),
2316    'mm2(100)' :(12,14, 0,-1),'mm2(010)' :(13,15, 0,-1),'mm2(110)' :( 6,11, 0,-1),'mm2(120)' :(15,14, 0,-1),
2317    'mm2(210)' :(16,15, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222(100)' :( 1,14, 0,-1),'222(010)' :( 1,15, 0,-1),
2318    '222(110)' :( 1,11, 0,-1),'2/m(001)' :( 1,16,14,-1),'2/m(100)' :( 1,25,12,-1),'2/m(010)' :( 1,28,13,-1),
2319    '2/m(110)' :( 1,19, 6,-1),'2/m(120)' :( 1,27,15,-1),'2/m(210)' :( 1,26,16,-1),'2/m(+-0)' :( 1,20,17,-1),
2320    'm(001)'   :(23,16,14,23),'m(100)'   :(26,25,12,26),'m(010)'   :(27,28,13,27),'m(110)'   :(18,19, 6,18),
2321    'm(120)'   :(24,27,15,24),'m(210)'   :(25,26,16,25),'m(+-0)'   :(17,20, 7,17),'2(001)'   :(14,16,14,23),
2322    '2(100)'   :(12,25,12,26),'2(010)'   :(13,28,13,27),'2(110)'   :( 6,19, 6,18),'2(120)'   :(15,27,15,24),
2323    '2(210)'   :(16,26,16,25),'2(+-0)'   :( 7,20, 7,17),'-1'       :( 1,29,28, 0)
2324    }
2325       
2326CSxinel = [[],      # 0th empty - indices are Fortran style
2327    [[0,0,0],[ 0.0, 0.0, 0.0]],      #1  0  0  0
2328    [[1,1,1],[ 1.0, 1.0, 1.0]],      #2  X  X  X
2329    [[1,1,1],[ 1.0, 1.0,-1.0]],      #3  X  X -X
2330    [[1,1,1],[ 1.0,-1.0, 1.0]],      #4  X -X  X
2331    [[1,1,1],[ 1.0,-1.0,-1.0]],      #5 -X  X  X
2332    [[1,1,0],[ 1.0, 1.0, 0.0]],      #6  X  X  0
2333    [[1,1,0],[ 1.0,-1.0, 0.0]],      #7  X -X  0
2334    [[1,0,1],[ 1.0, 0.0, 1.0]],      #8  X  0  X
2335    [[1,0,1],[ 1.0, 0.0,-1.0]],      #9  X  0 -X
2336    [[0,1,1],[ 0.0, 1.0, 1.0]],      #10  0  Y  Y
2337    [[0,1,1],[ 0.0, 1.0,-1.0]],      #11 0  Y -Y
2338    [[1,0,0],[ 1.0, 0.0, 0.0]],      #12  X  0  0
2339    [[0,1,0],[ 0.0, 1.0, 0.0]],      #13  0  Y  0
2340    [[0,0,1],[ 0.0, 0.0, 1.0]],      #14  0  0  Z
2341    [[1,1,0],[ 1.0, 2.0, 0.0]],      #15  X 2X  0
2342    [[1,1,0],[ 2.0, 1.0, 0.0]],      #16 2X  X  0
2343    [[1,1,2],[ 1.0, 1.0, 1.0]],      #17  X  X  Z
2344    [[1,1,2],[ 1.0,-1.0, 1.0]],      #18  X -X  Z
2345    [[1,2,1],[ 1.0, 1.0, 1.0]],      #19  X  Y  X
2346    [[1,2,1],[ 1.0, 1.0,-1.0]],      #20  X  Y -X
2347    [[1,2,2],[ 1.0, 1.0, 1.0]],      #21  X  Y  Y
2348    [[1,2,2],[ 1.0, 1.0,-1.0]],      #22  X  Y -Y
2349    [[1,2,0],[ 1.0, 1.0, 0.0]],      #23  X  Y  0
2350    [[1,0,2],[ 1.0, 0.0, 1.0]],      #24  X  0  Z
2351    [[0,1,2],[ 0.0, 1.0, 1.0]],      #25  0  Y  Z
2352    [[1,1,2],[ 1.0, 2.0, 1.0]],      #26  X 2X  Z
2353    [[1,1,2],[ 2.0, 1.0, 1.0]],      #27 2X  X  Z
2354    [[1,2,3],[ 1.0, 1.0, 1.0]],      #28  X  Y  Z
2355    ]
2356
2357CSuinel = [[],      # 0th empty - indices are Fortran style
2358    [[1,1,1,0,0,0],[ 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],[1,0,0,0,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #1  A  A  A  0  0  0
2359    [[1,1,2,0,0,0],[ 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],[1,0,1,0,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #2  A  A  C  0  0  0
2360    [[1,2,1,0,0,0],[ 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],[1,1,0,0,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #3  A  B  A  0  0  0
2361    [[1,2,2,0,0,0],[ 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],[1,1,0,0,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #4  A  B  B  0  0  0
2362    [[1,1,1,2,2,2],[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],[1,0,0,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #5  A  A  A  D  D  D
2363    [[1,1,1,2,2,2],[ 1.0, 1.0, 1.0, 1.0,-1.0,-1.0],[1,0,0,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #6  A  A  A  D -D -D
2364    [[1,1,1,2,2,2],[ 1.0, 1.0, 1.0, 1.0,-1.0, 1.0],[1,0,0,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #7  A  A  A  D -D  D
2365    [[1,1,1,2,2,2],[ 1.0, 1.0, 1.0, 1.0, 1.0,-1.0],[1,0,0,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #8  A  A  A  D  D -D
2366    [[1,1,2,1,0,0],[ 1.0, 1.0, 1.0, 0.5, 0.0, 0.0],[1,0,1,0,0,0],[1.0,1.0,1.0,0.5,0.0,0.0]],    #9  A  A  C A/2 0  0
2367    [[1,2,3,0,0,0],[ 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],[1,1,1,0,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #10  A  B  C  0  0  0
2368    [[1,1,2,3,0,0],[ 1.0, 1.0, 1.0, 1.0, 0.0, 0.0],[1,0,1,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #11  A  A  C  D  0  0
2369    [[1,2,1,0,3,0],[ 1.0, 1.0, 1.0, 0.0, 1.0, 0.0],[1,1,0,0,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #12  A  B  A  0  E  0
2370    [[1,2,2,0,0,3],[ 1.0, 1.0, 1.0, 0.0, 0.0, 1.0],[1,1,0,0,0,1],[1.0,1.0,1.0,0.0,0.0,0.0]],    #13  A  B  B  0  0  F
2371    [[1,2,3,2,0,0],[ 1.0, 1.0, 1.0, 0.5, 0.0, 0.0],[1,1,1,0,0,0],[1.0,1.0,1.0,0.0,0.5,0.0]],    #14  A  B  C B/2 0  0
2372    [[1,2,3,1,0,0],[ 1.0, 1.0, 1.0, 0.5, 0.0, 0.0],[1,1,1,0,0,0],[1.0,1.0,1.0,0.0,0.5,0.0]],    #15  A  B  C A/2 0  0
2373    [[1,2,3,4,0,0],[ 1.0, 1.0, 1.0, 1.0, 0.0, 0.0],[1,1,1,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #16  A  B  C  D  0  0
2374    [[1,2,3,0,4,0],[ 1.0, 1.0, 1.0, 0.0, 1.0, 0.0],[1,1,1,0,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #17  A  B  C  0  E  0
2375    [[1,2,3,0,0,4],[ 1.0, 1.0, 1.0, 0.0, 0.0, 1.0],[1,1,1,0,0,1],[1.0,1.0,1.0,0.0,0.0,0.0]],    #18  A  B  C  0  0  F
2376    [[1,1,2,3,4,4],[ 1.0, 1.0, 1.0, 1.0, 1.0,-1.0],[1,0,1,1,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #19  A  A  C  D  E -E
2377    [[1,1,2,3,4,4],[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],[1,0,1,1,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #20  A  A  C  D  E  E
2378    [[1,2,1,3,4,3],[ 1.0, 1.0, 1.0, 1.0, 1.0,-1.0],[1,1,0,1,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #21  A  B  A  D  E -D
2379    [[1,2,1,3,4,3],[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],[1,1,0,1,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #22  A  B  A  D  E  D
2380    [[1,2,2,3,3,4],[ 1.0, 1.0, 1.0, 1.0,-1.0, 1.0],[1,1,0,1,0,1],[1.0,1.0,1.0,0.0,0.0,0.0]],    #23  A  B  B  D -D  F
2381    [[1,2,2,3,3,4],[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],[1,1,0,1,0,1],[1.0,1.0,1.0,0.0,0.0,0.0]],    #24  A  B  B  D  D  F
2382    [[1,2,3,2,4,4],[ 1.0, 1.0, 1.0, 0.5, 1.0, 2.0],[1,1,1,0,0,1],[1.0,1.0,1.0,0.5,0.0,0.0]],    #25  A  B  C B/2 F/2 F
2383    [[1,2,3,1,0,4],[ 1.0, 1.0, 1.0, 0.5, 0.0, 1.0],[1,1,1,0,0,1],[1.0,1.0,1.0,0.5,0.0,0.0]],    #26  A  B  C A/2  0  F
2384    [[1,2,3,2,4,0],[ 1.0, 1.0, 1.0, 0.5, 1.0, 0.0],[1,1,1,0,1,0],[1.0,1.0,1.0,0.5,0.0,0.0]],    #27  A  B  C B/2  E  0
2385    [[1,2,3,1,4,4],[ 1.0, 1.0, 1.0, 0.5, 1.0, 0.5],[1,1,1,0,1,0],[1.0,1.0,1.0,0.5,0.0,0.0]],    #28  A  B  C A/2  E E/2
2386    [[1,2,3,4,5,6],[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],[1,1,1,1,1,1],[1.0,1.0,1.0,0.0,0.0,0.0]],    #29  A  B  C  D  E   F
2387    ]
2388   
2389################################################################################
2390#### Site symmetry routines
2391################################################################################
2392   
2393def GetOprPtrName(key):
2394    'Needs a doc string'
2395    try:
2396        oprName = OprName[key][0]
2397    except KeyError:
2398        return key
2399    return oprName.replace('(','').replace(')','')
2400
2401def GetOprPtrNumber(key):
2402    'Needs a doc string'
2403    try:
2404        return OprName[key][1]
2405    except KeyError:
2406        return key
2407
2408def GetOprName(key):
2409    'Needs a doc string'
2410    return OprName[key][0]
2411
2412def GetKNsym(key):
2413    'Needs a doc string'
2414    try:
2415        return KNsym[key].strip()
2416    except KeyError:
2417        return 'sp'
2418
2419def GetNXUPQsym(siteSym):
2420    '''       
2421    The codes XUPQ are for lookup of symmetry constraints for position(X), thermal parm(U) & magnetic moments (P & Q)
2422    '''
2423    return NXUPQsym[siteSym]
2424
2425def GetCSxinel(siteSym): 
2426    "returns Xyz terms, multipliers, GUI flags"
2427    indx = GetNXUPQsym(siteSym.strip())
2428    return CSxinel[indx[0]]
2429   
2430def GetCSuinel(siteSym):
2431    "returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
2432    indx = GetNXUPQsym(siteSym.strip())
2433    return CSuinel[indx[1]]
2434   
2435def GetCSpqinel(SpnFlp,dupDir): 
2436    "returns Mxyz terms, multipliers, GUI flags"
2437    CSI = [[1,2,3],[1.0,1.0,1.0]]
2438    for sopr in dupDir:
2439#        print (sopr,dupDir[sopr])
2440        opr = sopr.replace("'",'')
2441        indx = GetNXUPQsym(opr)
2442        if SpnFlp[dupDir[sopr]] > 0:
2443            csi = CSxinel[indx[2]]  #P
2444        else:
2445            csi = CSxinel[indx[3]]  #Q
2446#        print(opr,indx,csi,CSI)
2447        if not len(csi):
2448            return [[0,0,0],[0.,0.,0.]]
2449        for kcs in [0,1,2]:
2450            if csi[0][kcs] == 0 and CSI[0][kcs] != 0:
2451                jcs = CSI[0][kcs]
2452                for ics in [0,1,2]:
2453                    if CSI[0][ics] == jcs:
2454                        CSI[0][ics] = 0
2455                        CSI[1][ics] = 0.
2456                    elif CSI[0][ics] > jcs:
2457                        CSI[0][ics] = CSI[0][ics]-1
2458            elif (CSI[0][kcs] == csi[0][kcs]) and (CSI[1][kcs] != csi[1][kcs]):
2459                CSI[1][kcs] = csi[1][kcs]
2460            elif CSI[0][kcs] >= csi[0][kcs]:
2461                CSI[0][kcs] = min(CSI[0][kcs],csi[0][kcs])
2462                if CSI[1][kcs] != csi[1][kcs]:
2463                    if CSI[1][kcs] == 1.:
2464                        CSI[1][kcs] = csi[1][kcs]
2465#        print(CSI)
2466    return CSI
2467   
2468def getTauT(tau,sop,ssop,XYZ,wave=np.zeros(3)):
2469    phase = np.sum(XYZ*wave)
2470    ssopinv = nl.inv(ssop[0])
2471    mst = ssopinv[3][:3]
2472    epsinv = ssopinv[3][3]
2473    sdet = nl.det(sop[0])
2474    ssdet = nl.det(ssop[0])
2475    dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
2476    dT = 1.0
2477    if np.any(dtau%.5):
2478        sumdtau = np.sum(dtau%.5)
2479        dT = 0.
2480        if np.abs(sumdtau-.5) > 1.e-4:
2481            dT = np.tan(np.pi*sumdtau)
2482    tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]+phase)
2483    return sdet,ssdet,dtau,dT,tauT
2484   
2485def OpsfromStringOps(A,SGData,SSGData):
2486    SGOps = SGData['SGOps']
2487    SSGOps = SSGData['SSGOps']
2488    Ax = A.split('+')
2489    Ax[0] = int(Ax[0])
2490    iC = 1
2491    if Ax[0] < 0:
2492        iC = -1
2493    Ax[0] = abs(Ax[0])
2494    nA = Ax[0]%100-1
2495    nC = Ax[0]//100
2496    unit = [0,0,0]
2497    if len(Ax) > 1:
2498        unit = eval('['+Ax[1]+']')
2499    return SGOps[nA],SSGOps[nA],iC,SGData['SGCen'][nC],unit
2500   
2501def GetSSfxuinel(waveType,Stype,nH,XYZ,SGData,SSGData,debug=False):
2502   
2503    def orderParms(CSI):
2504        parms = [0,]
2505        for csi in CSI:
2506            for i in [0,1,2]:
2507                if csi[i] not in parms:
2508                    parms.append(csi[i])
2509        for csi in CSI:
2510            for i in [0,1,2]:
2511                csi[i] = parms.index(csi[i])
2512        return CSI
2513       
2514    def fracCrenel(tau,Toff,Twid):
2515        Tau = (tau-Toff[:,nxs])%1.
2516        A = np.where(Tau<Twid[:,nxs],1.,0.)
2517        return A
2518       
2519    def fracFourier(tau,nH,fsin,fcos):
2520        SA = np.sin(2.*nH*np.pi*tau)
2521        CB = np.cos(2.*nH*np.pi*tau)
2522        A = SA[nxs,nxs,:]*fsin[:,:,nxs]
2523        B = CB[nxs,nxs,:]*fcos[:,:,nxs]
2524        return A+B
2525       
2526    def posFourier(tau,nH,psin,pcos):
2527        SA = np.sin(2*nH*np.pi*tau)
2528        CB = np.cos(2*nH*np.pi*tau)
2529        A = SA[nxs,nxs,:]*psin[:,:,nxs]
2530        B = CB[nxs,nxs,:]*pcos[:,:,nxs]
2531        return A+B   
2532
2533    def posSawtooth(tau,Toff,slopes):
2534        Tau = (tau-Toff)%1.
2535        A = slopes[:,nxs]*Tau
2536        return A
2537   
2538    def posZigZag(tau,Tmm,XYZmax):
2539        DT = Tmm[1]-Tmm[0]
2540        slopeUp = 2.*XYZmax/DT
2541        slopeDn = 2.*XYZmax/(1.-DT)
2542        A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax+slopeUp*((t-Tmm[0])%1.),XYZmax-slopeDn*((t-Tmm[1])%1.)) for t in tau])
2543        return A
2544
2545    def posBlock(tau,Tmm,XYZmax):
2546        A = np.array([np.where(Tmm[0] < t <= Tmm[1],XYZmax,-XYZmax) for t in tau])
2547        return A
2548       
2549    def DoFrac():
2550        dF = None
2551        dFTP = None
2552        if siteSym == '1':
2553            CSI = [[1,0],[2,0]],2*[[1.,0.],]
2554        elif siteSym == '-1':
2555            CSI = [[1,0],[0,0]],2*[[1.,0.],]
2556        else:
2557            delt2 = np.eye(2)*0.001
2558            FSC = np.ones(2,dtype='i')
2559            CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
2560            if 'Crenel' in waveType:
2561                dF = np.zeros_like(tau)
2562            else:
2563                dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2564            dFT = np.zeros_like(dF)
2565            dFTP = []
2566            for i in SdIndx:
2567                sop = Sop[i]
2568                ssop = SSop[i]           
2569                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2570                fsc = np.ones(2,dtype='i')
2571                if 'Crenel' in waveType:
2572                    dFT = np.zeros_like(tau)
2573                    fsc = [1,1]
2574                else:   #Fourier
2575                    dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
2576                    dFT = nl.det(sop[0])*dFT
2577                    dFT = dFT[:,np.argsort(tauT)]
2578                    dFT[0] *= ssdet
2579                    dFT[1] *= sdet
2580                    dFTP.append(dFT)
2581               
2582                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2583                        fsc = [1,1]
2584                        if dT:
2585                            CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
2586                        else:
2587                            CSI = [[[1,0],[0,0]],[[1.,0.],[0.,0.]]]
2588                        FSC = np.zeros(2,dtype='i')
2589                        return CSI,dF,dFTP
2590                    else:
2591                        for i in range(2):
2592                            if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
2593                                fsc[i] = 1
2594                            else:
2595                                fsc[i] = 0
2596                        FSC &= fsc
2597                        if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc)
2598            n = -1
2599            for i,F in enumerate(FSC):
2600                if F:
2601                    n += 1
2602                    CSI[0][i] = n+1
2603                    CSI[1][i] = 1.0
2604           
2605        return CSI,dF,dFTP
2606       
2607    def DoXYZ():
2608        dX,dXTP = None,None
2609        if siteSym == '1':
2610            CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2611        elif siteSym == '-1':
2612            CSI = [[1,0,0],[2,0,0],[3,0,0], [0,0,0],[0,0,0],[0,0,0]],3*[[1.,0.,0.],]+3*[[0.,0.,0.],]
2613        else:
2614            delt4 = np.ones(4)*0.001
2615            delt5 = np.ones(5)*0.001
2616            delt6 = np.eye(6)*0.001
2617            if 'Fourier' in waveType:
2618                dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2619                  #3x6x12 modulated position array (X,Spos,tau)& force positive
2620                CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2621            elif waveType == 'Sawtooth':
2622                dX = posSawtooth(tau,delt4[0],delt4[1:])
2623                CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0]]),
2624                    np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2625            elif waveType in ['ZigZag','Block']:
2626                if waveType == 'ZigZag':
2627                    dX = posZigZag(tau,delt5[:2],delt5[2:])
2628                else:
2629                    dX = posBlock(tau,delt5[:2],delt5[2:])
2630                CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2631                    np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2632            XSC = np.ones(6,dtype='i')
2633            dXTP = []
2634            for i in SdIndx:
2635                sop = Sop[i]
2636                ssop = SSop[i]
2637                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2638                xsc = np.ones(6,dtype='i')
2639                if 'Fourier' in waveType:
2640                    dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,nxs,nxs]
2641                elif waveType == 'Sawtooth':
2642                    dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,nxs,nxs]
2643                elif waveType == 'ZigZag':
2644                    dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2645                elif waveType == 'Block':
2646                    dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2647                dXT = np.inner(sop[0],dXT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2648                dXT = np.swapaxes(dXT,1,2)      # back to array(3x6x49)
2649                dXT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2650                dXTP.append(dXT)
2651                if waveType == 'Fourier':
2652                    for i in range(3):
2653                        if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2654                            xsc[i] = 0
2655                        if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2656                            xsc[i+3] = 0
2657                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2658                        xsc[3:6] = 0
2659                        CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2660                            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2661                        if dT:
2662                            if '(x)' in siteSym:
2663                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2664                                if 'm' in siteSym and len(SdIndx) == 1:
2665                                    CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2666                            elif '(y)' in siteSym:
2667                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2668                                if 'm' in siteSym and len(SdIndx) == 1:
2669                                    CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2670                            elif '(z)' in siteSym:
2671                                CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2672                                if 'm' in siteSym and len(SdIndx) == 1:
2673                                    CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2674                        else:
2675                            CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2676                    if '4/mmm' in laue:
2677                        if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2678                            if '(xy)' in siteSym:
2679                                CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2680                                if dT:
2681                                    CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2682                                else:
2683                                    CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2684                        if '(xy)' in siteSym or '(+-0)' in siteSym:
2685                            mul = 1
2686                            if '(+-0)' in siteSym:
2687                                mul = -1
2688                            if np.allclose(dX[0,0,:],dXT[1,0,:]):
2689                                CSI[0][3:5] = [[11,0,0],[11,0,0]]
2690                                CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2691                                xsc[3:5] = 0
2692                            if np.allclose(dX[0,3,:],dXT[0,4,:]):
2693                                CSI[0][:2] = [[12,0,0],[12,0,0]]
2694                                CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2695                                xsc[:2] = 0
2696                XSC &= xsc
2697                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2698            if waveType == 'Fourier':
2699                n = -1
2700                if debug: print (XSC)
2701                for i,X in enumerate(XSC):
2702                    if X:
2703                        n += 1
2704                        CSI[0][i][0] = n+1
2705                        CSI[1][i][0] = 1.0
2706           
2707        return list(CSI),dX,dXTP
2708       
2709    def DoUij():
2710        dU,dUTP = None,None
2711        if siteSym == '1':
2712            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2713                [7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2714        elif siteSym == '-1':
2715            CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],   \
2716                6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2717        else:
2718            tau = np.linspace(0,1,49,True)
2719            delt12 = np.eye(12)*0.0001
2720            dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
2721            CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2722            USC = np.ones(12,dtype='i')
2723            dUTP = []
2724            dtau = 0.
2725            for i in SdIndx:
2726                sop = Sop[i]
2727                ssop = SSop[i]
2728                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2729                usc = np.ones(12,dtype='i')
2730                dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
2731                dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3)    #convert dUT to 12x49x3x3
2732                dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2733                dUT = np.array(U2Uij(dUijT))    #convert to 6x12x49
2734                dUT = dUT[:,:,np.argsort(tauT)]
2735                dUT[:,:6,:] *=(ssdet*sdet)
2736                dUTP.append(dUT)
2737                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2738                    if dT:
2739                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2740                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2741                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2742                        [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2743                    else:
2744                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2745                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2746                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2747                        [0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2748                    if 'mm2(x)' in siteSym and dT:
2749                        CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2750                        USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2751                    elif '(xy)' in siteSym and dT:
2752                        CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2753                            [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2754                        CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2755                        USC = [1,1,1,1,1,1,1,1,1,1,1,1]                             
2756                    elif '(x)' in siteSym and dT:
2757                        CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2758                    elif '(y)' in siteSym and dT:
2759                        CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2760                    elif '(z)' in siteSym and dT:
2761                        CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2762                    for i in range(6):
2763                        if not USC[i]:
2764                            CSI[0][i] = [0,0,0]
2765                            CSI[1][i] = [0.,0.,0.]
2766                            CSI[0][i+6] = [0,0,0]
2767                            CSI[1][i+6] = [0.,0.,0.]
2768                else:                       
2769                    for i in range(6):
2770                        if not np.allclose(dU[i,i,:],dUT[i,i,:]):  #sin part
2771                            usc[i] = 0
2772                        if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]):   #cos part
2773                            usc[i+6] = 0
2774                    if np.any(dUT[1,0,:]):
2775                        if '4/m' in siteSym:
2776                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2777                            if ssop[1][3]:
2778                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2779                                usc[9] = 1
2780                            else:
2781                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2782                                usc[9] = 0
2783                        elif '4' in siteSym:
2784                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2785                            CSI[0][:2] = [[11,0,0],[11,0,0]]
2786                            if ssop[1][3]:
2787                                CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2788                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2789                                usc[2] = 0
2790                                usc[8] = 0
2791                                usc[3] = 1
2792                                usc[9] = 1
2793                            else:
2794                                CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2795                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2796                                usc[2] = 1
2797                                usc[8] = 1
2798                                usc[3] = 0               
2799                                usc[9] = 0
2800                        elif 'xy' in siteSym or '+-0' in siteSym:
2801                            if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2802                                CSI[0][4:6] = [[12,0,0],[12,0,0]]
2803                                CSI[0][6:8] = [[11,0,0],[11,0,0]]
2804                                CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2805                                CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2806                                usc[4:6] = 0
2807                                usc[6:8] = 0
2808                           
2809                    if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2810                USC &= usc
2811            if debug: print (USC)
2812            if not np.any(dtau%.5):
2813                n = -1
2814                for i,U in enumerate(USC):
2815                    if U:
2816                        n += 1
2817                        CSI[0][i][0] = n+1
2818                        CSI[1][i][0] = 1.0
2819   
2820        return list(CSI),dU,dUTP
2821   
2822    def DoMag():
2823        CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2824        if siteSym == '1':
2825            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2826        elif siteSym in ['-1','mmm','2/m(x)','2/m(y)','2/m(z)','4/mmm001']:
2827            CSI = 3*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3*[[0.,0.,0.],]+3*[[1.,0.,0.],]
2828        else:
2829            delt6 = np.eye(6)*0.001
2830            dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,nxs,nxs]
2831              #3x6x12 modulated moment array (M,Spos,tau)& force positive
2832            CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2833            MSC = np.ones(6,dtype='i')
2834            dMTP = []
2835            for i in SdIndx:
2836                sop = Sop[i]
2837                ssop = SSop[i]
2838                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2839                msc = np.ones(6,dtype='i')
2840                dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,nxs,nxs]
2841                dMT = np.inner(sop[0],dMT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2842                dMT = np.swapaxes(dMT,1,2)      # back to array(3x6x49)
2843                dMT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2844                dMTP.append(dMT)
2845                for i in range(3):
2846                    if not np.allclose(dM[i,i,:],dMT[i,i,:]):
2847                        msc[i] = 0
2848                    if not np.allclose(dM[i,i+3,:],dMT[i,i+3,:]):
2849                        msc[i+3] = 0
2850                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2851                    msc[3:6] = 0
2852                    CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2853                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2854                    if dT:
2855                        if '(x)' in siteSym:
2856                            CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2857                            if 'm' in siteSym and len(SdIndx) == 1:
2858                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2859                        elif '(y)' in siteSym:
2860                            CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2861                            if 'm' in siteSym and len(SdIndx) == 1:
2862                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2863                        elif '(z)' in siteSym:
2864                            CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2865                            if 'm' in siteSym and len(SdIndx) == 1:
2866                                CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2867                    else:
2868                        CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2869                if '4/mmm' in laue:
2870                    if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2871                        if '(xy)' in siteSym:
2872                            CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2873                            if dT:
2874                                CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2875                            else:
2876                                CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2877                    if '(xy)' in siteSym or '(+-0)' in siteSym:
2878                        mul = 1
2879                        if '(+-0)' in siteSym:
2880                            mul = -1
2881                        if np.allclose(dM[0,0,:],dMT[1,0,:]):
2882                            CSI[0][3:5] = [[11,0,0],[11,0,0]]
2883                            CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2884                            msc[3:5] = 0
2885                        if np.allclose(dM[0,3,:],dMT[0,4,:]):
2886                            CSI[0][:2] = [[12,0,0],[12,0,0]]
2887                            CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2888                            msc[:2] = 0
2889                MSC &= msc
2890                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2891            n = -1
2892            if debug: print (MSC)
2893            for i,M in enumerate(MSC):
2894                if M:
2895                    n += 1
2896                    CSI[0][i][0] = n+1
2897                    CSI[1][i][0] = 1.0
2898
2899        return CSI,None,None
2900       
2901    if debug: print ('super space group: '+SSGData['SSpGrp'])
2902    xyz = np.array(XYZ)%1.
2903    SGOps = copy.deepcopy(SGData['SGOps'])
2904    laue = SGData['SGLaue']
2905    siteSym = SytSym(XYZ,SGData)[0].strip()
2906    if debug: print ('siteSym: '+siteSym)
2907    SSGOps = copy.deepcopy(SSGData['SSGOps'])
2908    #expand ops to include inversions if any
2909    if SGData['SGInv']:
2910        for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
2911            SGOps.append([-op[0],-op[1]%1.])
2912            SSGOps.append([-sop[0],-sop[1]%1.])
2913    #build set of sym ops around special position       
2914    SSop = []
2915    Sop = []
2916    Sdtau = []
2917    for iop,Op in enumerate(SGOps):         
2918        nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
2919        if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
2920            SSop.append(SSGOps[iop])
2921            Sop.append(SGOps[iop])
2922            ssopinv = nl.inv(SSGOps[iop][0])
2923            mst = ssopinv[3][:3]
2924            epsinv = ssopinv[3][3]
2925            Sdtau.append(np.sum(mst*(XYZ-SGOps[iop][1])-epsinv*SSGOps[iop][1][3]))
2926    SdIndx = np.argsort(np.array(Sdtau))     # just to do in sensible order
2927    if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
2928    #setup displacement arrays
2929    tau = np.linspace(-1,1,49,True)
2930    #make modulation arrays - one parameter at a time
2931    if Stype == 'Sfrac':
2932        CSI,dF,dFTP = DoFrac()
2933    elif Stype == 'Spos':
2934        CSI,dF,dFTP = DoXYZ()
2935        CSI[0] = orderParms(CSI[0])
2936    elif Stype == 'Sadp':
2937        CSI,dF,dFTP = DoUij()
2938        CSI[0] = orderParms(CSI[0]) 
2939    elif Stype == 'Smag':
2940        CSI,dF,dFTP = DoMag()
2941       
2942    if debug:
2943        return CSI,dF,dFTP
2944    else:
2945        return CSI
2946   
2947def MustrainNames(SGData):
2948    'Needs a doc string'
2949    laue = SGData['SGLaue']
2950    uniq = SGData['SGUniq']
2951    if laue in ['m3','m3m']:
2952        return ['S400','S220']
2953    elif laue in ['6/m','6/mmm','3m1']:
2954        return ['S400','S004','S202']
2955    elif laue in ['31m','3']:
2956        return ['S400','S004','S202','S211']
2957    elif laue in ['3R','3mR']:
2958        return ['S400','S220','S310','S211']
2959    elif laue in ['4/m','4/mmm']:
2960        return ['S400','S004','S220','S022']
2961    elif laue in ['mmm']:
2962        return ['S400','S040','S004','S220','S202','S022']
2963    elif laue in ['2/m']:
2964        SHKL = ['S400','S040','S004','S220','S202','S022']
2965        if uniq == 'a':
2966            SHKL += ['S013','S031','S211']
2967        elif uniq == 'b':
2968            SHKL += ['S301','S103','S121']
2969        elif uniq == 'c':
2970            SHKL += ['S130','S310','S112']
2971        return SHKL
2972    else:
2973        SHKL = ['S400','S040','S004','S220','S202','S022']
2974        SHKL += ['S310','S103','S031','S130','S301','S013']
2975        SHKL += ['S211','S121','S112']
2976        return SHKL
2977       
2978def HStrainVals(HSvals,SGData):
2979    laue = SGData['SGLaue']
2980    uniq = SGData['SGUniq']
2981    DIJ = np.zeros(6)
2982    if laue in ['m3','m3m']:
2983        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
2984    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2985        DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
2986    elif laue in ['3R','3mR']:
2987        DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
2988    elif laue in ['4/m','4/mmm']:
2989        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
2990    elif laue in ['mmm']:
2991        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2992    elif laue in ['2/m']:
2993        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2994        if uniq == 'a':
2995            DIJ[5] = HSvals[3]
2996        elif uniq == 'b':
2997            DIJ[4] = HSvals[3]
2998        elif uniq == 'c':
2999            DIJ[3] = HSvals[3]
3000    else:
3001        DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
3002    return DIJ
3003
3004def HStrainNames(SGData):
3005    'Needs a doc string'
3006    laue = SGData['SGLaue']
3007    uniq = SGData['SGUniq']
3008    if laue in ['m3','m3m']:
3009        return ['D11','eA']         #add cubic strain term
3010    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3011        return ['D11','D33']
3012    elif laue in ['3R','3mR']:
3013        return ['D11','D12']
3014    elif laue in ['4/m','4/mmm']:
3015        return ['D11','D33']
3016    elif laue in ['mmm']:
3017        return ['D11','D22','D33']
3018    elif laue in ['2/m']:
3019        Dij = ['D11','D22','D33']
3020        if uniq == 'a':
3021            Dij += ['D23']
3022        elif uniq == 'b':
3023            Dij += ['D13']
3024        elif uniq == 'c':
3025            Dij += ['D12']
3026        return Dij
3027    else:
3028        Dij = ['D11','D22','D33','D12','D13','D23']
3029        return Dij
3030   
3031def MustrainCoeff(HKL,SGData):
3032    'Needs a doc string'
3033    #NB: order of terms is the same as returned by MustrainNames
3034    laue = SGData['SGLaue']
3035    uniq = SGData['SGUniq']
3036    h,k,l = HKL
3037    Strm = []
3038    if laue in ['m3','m3m']:
3039        Strm.append(h**4+k**4+l**4)
3040        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
3041    elif laue in ['6/m','6/mmm','3m1']:
3042        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
3043        Strm.append(l**4)
3044        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
3045    elif laue in ['31m','3']:
3046        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
3047        Strm.append(l**4)
3048        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
3049        Strm.append(4.0*h*k*l*(h+k))
3050    elif laue in ['3R','3mR']:
3051        Strm.append(h**4+k**4+l**4)
3052        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
3053        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
3054        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
3055    elif laue in ['4/m','4/mmm']:
3056        Strm.append(h**4+k**4)
3057        Strm.append(l**4)
3058        Strm.append(3.0*(h*k)**2)
3059        Strm.append(3.0*((h*l)**2+(k*l)**2))
3060    elif laue in ['mmm']:
3061        Strm.append(h**4)
3062        Strm.append(k**4)
3063        Strm.append(l**4)
3064        Strm.append(3.0*(h*k)**2)
3065        Strm.append(3.0*(h*l)**2)
3066        Strm.append(3.0*(k*l)**2)
3067    elif laue in ['2/m']:
3068        Strm.append(h**4)
3069        Strm.append(k**4)
3070        Strm.append(l**4)
3071        Strm.append(3.0*(h*k)**2)
3072        Strm.append(3.0*(h*l)**2)
3073        Strm.append(3.0*(k*l)**2)
3074        if uniq == 'a':
3075            Strm.append(2.0*k*l**3)
3076            Strm.append(2.0*l*k**3)
3077            Strm.append(4.0*k*l*h**2)
3078        elif uniq == 'b':
3079            Strm.append(2.0*l*h**3)
3080            Strm.append(2.0*h*l**3)
3081            Strm.append(4.0*h*l*k**2)
3082        elif uniq == 'c':
3083            Strm.append(2.0*h*k**3)
3084            Strm.append(2.0*k*h**3)
3085            Strm.append(4.0*h*k*l**2)
3086    else:
3087        Strm.append(h**4)
3088        Strm.append(k**4)
3089        Strm.append(l**4)
3090        Strm.append(3.0*(h*k)**2)
3091        Strm.append(3.0*(h*l)**2)
3092        Strm.append(3.0*(k*l)**2)
3093        Strm.append(2.0*k*h**3)
3094        Strm.append(2.0*h*l**3)
3095        Strm.append(2.0*l*k**3)
3096        Strm.append(2.0*h*k**3)
3097        Strm.append(2.0*l*h**3)
3098        Strm.append(2.0*k*l**3)
3099        Strm.append(4.0*k*l*h**2)
3100        Strm.append(4.0*h*l*k**2)
3101        Strm.append(4.0*k*h*l**2)
3102    return Strm
3103
3104def MuShklMean(SGData,Amat,Shkl):
3105   
3106    def genMustrain(xyz,Shkl):
3107        uvw = np.inner(Amat.T,xyz)
3108        Strm = np.array(MustrainCoeff(uvw,SGData))
3109        Sum = np.sum(np.multiply(Shkl,Strm))
3110        Sum = np.where(Sum > 0.01,Sum,0.01)
3111        Sum = np.sqrt(Sum)
3112        return Sum*xyz
3113       
3114    PHI = np.linspace(0.,360.,30,True)
3115    PSI = np.linspace(0.,180.,30,True)
3116    X = np.outer(npcosd(PHI),npsind(PSI))
3117    Y = np.outer(npsind(PHI),npsind(PSI))
3118    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3119    XYZ = np.dstack((X,Y,Z))
3120    XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
3121    return np.sqrt(np.sum(XYZ**2)/900.)
3122   
3123def Muiso2Shkl(muiso,SGData,cell):
3124    "this is to convert isotropic mustrain to generalized Shkls"
3125    import GSASIIlattice as G2lat
3126    A = G2lat.cell2AB(cell)[0]
3127   
3128    def minMus(Shkl,muiso,H,SGData,A):
3129        U = np.inner(A.T,H)
3130        S = np.array(MustrainCoeff(U,SGData))
3131        Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,nxs]),axis=0))
3132        rad = np.sqrt(np.sum((Sum[:,nxs]*H)**2,axis=1))
3133        return (muiso-rad)**2
3134       
3135    laue = SGData['SGLaue']
3136    PHI = np.linspace(0.,360.,60,True)
3137    PSI = np.linspace(0.,180.,60,True)
3138    X = np.outer(npsind(PHI),npsind(PSI))
3139    Y = np.outer(npcosd(PHI),npsind(PSI))
3140    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3141    HKL = np.dstack((X,Y,Z))
3142    if laue in ['m3','m3m']:
3143        S0 = [1000.,1000.]
3144    elif laue in ['6/m','6/mmm','3m1']:
3145        S0 = [1000.,1000.,1000.]
3146    elif laue in ['31m','3']:
3147        S0 = [1000.,1000.,1000.,1000.]
3148    elif laue in ['3R','3mR']:
3149        S0 = [1000.,1000.,1000.,1000.]
3150    elif laue in ['4/m','4/mmm']:
3151        S0 = [1000.,1000.,1000.,1000.]
3152    elif laue in ['mmm']:
3153        S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
3154    elif laue in ['2/m']:
3155        S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
3156    else:
3157        S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000., 
3158            1000.,1000.,0.,0.,0.]
3159    S0 = np.array(S0)
3160    HKL = np.reshape(HKL,(-1,3))
3161    result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
3162    return result[0]
3163       
3164def PackRot(SGOps):
3165    IRT = []
3166    for ops in SGOps:
3167        M = ops[0]
3168        irt = 0
3169        for j in range(2,-1,-1):
3170            for k in range(2,-1,-1):
3171                irt *= 3
3172                irt += M[k][j]
3173        IRT.append(int(irt))
3174    return IRT
3175       
3176def SytSym(XYZ,SGData):
3177    '''
3178    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
3179
3180    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
3181    :param SGData: from SpcGroup
3182    :Returns: a four element tuple:
3183
3184     * The 1st element is a code for the site symmetry (see GetKNsym)
3185     * The 2nd element is the site multiplicity
3186     * Ndup number of overlapping operators
3187     * dupDir Dict - dictionary of overlapping operators
3188
3189    '''
3190    Mult = 1
3191    Isym = 0
3192    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
3193        Isym = 1073741824
3194    Jdup = 0
3195    Ndup = 0
3196    dupDir = {}
3197    inv = SGData['SGInv']+1
3198    icen = SGData['SGCen']
3199    if SGData['SGFixed']:       #already in list of operators
3200        inv = 1
3201    Xeqv = list(GenAtom(XYZ,SGData,True))
3202#    for xeqv in Xeqv:   print(xeqv)
3203    IRT = PackRot(SGData['SGOps'])
3204    L = -1
3205    for ic,cen in enumerate(icen):
3206        for invers in range(int(inv)):
3207            for io,ops in enumerate(SGData['SGOps']):
3208                irtx = (1-2*invers)*IRT[io]
3209                L += 1
3210                if not Xeqv[L][1]:
3211                    Ndup = io
3212                    Jdup += 1
3213                    jx = GetOprPtrNumber(str(irtx))   #[KN table no,op name,KNsym ptr]
3214                    if jx < 39:
3215                        px = GetOprName(str(irtx))
3216                        if Xeqv[L][-1] < 0:
3217                            if '(' in px:
3218                                px = px.split('(')
3219                                px[0] += "'"
3220                                px = '('.join(px)
3221                            else:   
3222                                px += "'"
3223                        dupDir[px] = L
3224                        Isym += 2**(jx-1)
3225    if Isym == 1073741824: Isym = 0
3226    Mult = len(SGData['SGOps'])*len(SGData['SGCen'])*inv//Jdup
3227         
3228    return GetKNsym(str(Isym)),Mult,Ndup,dupDir
3229   
3230def MagSytSym(SytSym,dupDir,SGData):
3231    '''
3232    site sym operations: 1,-1,2,3,-3,4,-4,6,-6,m need to be marked if spin inversion
3233    '''
3234    SGData['GenSym'],SGData['GenFlg'] = GetGenSym(SGData)[:2]
3235#    print('SGPtGrp',SGData['SGPtGrp'],'SytSym',SytSym,'MagSpGrp',SGData['MagSpGrp'])
3236#    print('dupDir',dupDir)
3237    SplitSytSym = SytSym.split('(')
3238    if SGData['SGGray']:
3239        return SytSym+"1'"
3240    if SytSym == '1':       #genersl position
3241        return SytSym
3242    if SplitSytSym[0] == SGData['SGPtGrp']:     #simple cases
3243        try:
3244            MagSytSym = SGData['MagSpGrp'].split()[1]
3245        except IndexError:
3246            MagSytSym = SGData['MagSpGrp'][1:].strip("1'")
3247        if len(SplitSytSym) > 1:
3248            MagSytSym += '('+SplitSytSym[1]
3249        return MagSytSym
3250    if len(dupDir) == 1:
3251        return list(dupDir.keys())[0]
3252   
3253   
3254    if '2/m' in SytSym:         #done I think; last 2wo might be not needed
3255        ops = {'(x)':['2(x)','m(x)'],'(y)':['2(y)','m(y)'],'(z)':['2(z)','m(z)'],
3256               '(100)':['2(100)','m(100)'],'(010)':['2(010)','m(010)'],'(001)':['2(001)','m(001)'],
3257               '(120)':['2(120)','m(120)'],'(210)':['2(210)','m(210)'],'(+-0)':['2(+-0)','m(+-0)'],
3258               '(110)':['2(110)','m(110)']}
3259   
3260    elif '4/mmm' in SytSym:
3261        ops = {'(x)':['4(x)','m(x)','m(y)','m(0+-)'],   #m(0+-) for cubic m3m?
3262               '(y)':['4(y)','m(y)','m(z)','m(+0-)'],   #m(+0-)
3263               '(z)':['4(z)','m(z)','m(x)','m(+-0)']}   #m(+-0)
3264    elif '4mm' in SytSym:
3265        ops = {'(x)':['4(x)','m(y)','m(yz)'],'(y)':['4(y)','m(z)','m(xz)'],'(z)':['4(z)','m(x)','m(110)']}
3266    elif '422' in SytSym:
3267        ops = {'(x)':['4(x)','2(y)','2(yz)'],'(y)':['4(y)','2(z)','2(xz)'],'(z)':['4(z)','2(x)','2(110)']}
3268    elif '-4m2' in SytSym:
3269        ops = {'(x)':['-4(x)','m(x)','2(yz)'],'(y)':['-4(y)','m(y)','2(xz)'],'(z)':['-4(z)','m(z)','2(110)']}
3270    elif '-42m' in SytSym:
3271        ops = {'(x)':['-4(x)','2(y)','m(yz)'],'(y)':['-4(y)','2(z)','m(xz)'],'(z)':['-4(z)','2(x)','m(110)']}
3272    elif '-4' in SytSym:
3273        ops = {'(x)':['-4(x)',],'(y)':['-4(y)',],'(z)':['-4(z)',],}
3274    elif '4' in SytSym:
3275        ops = {'(x)':['4(x)',],'(y)':['4(y)',],'(z)':['4(z)',],}
3276
3277    elif '222' in SytSym:
3278        ops = {'':['2(x)','2(y)','2(z)'],
3279                   '(x)':['2(y)','2(z)','2(x)'],'(y)':['2(x)','2(z)','2(y)'],'(z)':['2(x)','2(y)','2(z)'],
3280                   '(100)':['2(z)','2(100)','2(120)',],'(010)':['2(z)','2(010)','2(210)',],
3281                   '(110)':['2(z)','2(110)','2(+-0)',],}
3282    elif 'mm2' in SytSym:
3283        ops = {'(x)':['m(y)','m(z)','2(x)'],'(y)':['m(x)','m(z)','2(y)'],'(z)':['m(x)','m(y)','2(z)'],
3284               '(xy)':['m(+-0)','m(z)','2(110)'],'(yz)':['m(0+-)','m(xz)','2(yz)'],     #not 2(xy)!
3285               '(xz)':['m(+0-)','m(y)','2(xz)'],'(z100)':['m(100)','m(120)','2(z)'],
3286               '(z010)':['m(010)','m(210)','2(z)'],'(z110)':['m(110)','m(+-0)','2(z)'],
3287               '(+-0)':[ 'm(110)','m(z)','2(+-0)'],'(d100)':['m(yz)','m(0+-)','2(xz)'],
3288               '(d010)':['m(xz)','m(+0-)','2(y)'],'(d001)':['m(110)','m(+-0)','2(z)'],
3289               '(210)':['m(z)','m(010)','2(210)'],
3290               '(100)':['m(z)','m(120)','2(100)',],'(010)':['m(z)','m(210)','2(010)',],
3291               '(110)':['m(z)','m(+-0)','2(110)',],}
3292    elif 'mmm' in SytSym:
3293        ops = {'':['m(x)','m(y)','m(z)'],
3294                   '(100)':['m(z)','m(100)','m(120)',],'(010)':['m(z)','m(010)','m(210)',],
3295                   '(110)':['m(z)','m(110)','m(+-0)',],
3296                   '(x)':['m(x)','m(y)','m(z)'],'(y)':['m(x)','m(y)','m(z)'],'(z)':['m(x)','m(y)','m(z)'],}
3297       
3298    elif '32' in SytSym:
3299        ops = {'(120)':['3','2(120)',],'(100)':['3','2(100)']}
3300    elif '23' in SytSym:
3301        ops = {'':['2(x)','3(111)']}
3302    elif 'm3' in SytSym:
3303        ops = {'(100)':['(+-0)',],'(+--)':[],'(-+-)':[],'(--+)':[]}
3304    elif '3m' in SytSym:
3305        ops = {'(111)':['3(111)','m(+-0)',],'(+--)':['3(+--)','m(0+-)',],
3306               '(-+-)':['3(-+-)','m(+0-)',],'(--+)':['3(--+)','m(+-0)',],
3307               '(100)':['3','m(100)'],'(120)':['3','m(210)',]}
3308   
3309    if SytSym.split('(')[0] in ['6/m','6mm','-6m2','622','-6','-3','-3m',]:     #not simple cases
3310        MagSytSym = SytSym
3311        if "-1'" in dupDir:
3312            if '-6' in SytSym:
3313                MagSytSym = MagSytSym.replace('-6',"-6'")
3314            elif '-3m' in SytSym:
3315                MagSytSym = MagSytSym.replace('-3m',"-3'm'")
3316            elif '-3' in SytSym:
3317                MagSytSym = MagSytSym.replace('-3',"-3'")
3318        elif '-6m2' in SytSym:
3319            if "m'(110)" in dupDir:
3320                MagSytSym = "-6m'2'("+SytSym.split('(')[1]
3321        elif '6/m' in SytSym:
3322            if "m'(z)" in dupDir:
3323                MagSytSym = "6'/m'"
3324        elif '6mm' in SytSym:
3325            if "m'(110)" in dupDir:
3326                MagSytSym = "6'm'm"
3327        return MagSytSym
3328    try:
3329        axis = '('+SytSym.split('(')[1]
3330    except IndexError:
3331        axis = ''
3332    MagSytSym = ''
3333    for m in ops[axis]:
3334        if m in dupDir:
3335            MagSytSym += m.split('(')[0]
3336        else:
3337            MagSytSym += m.split('(')[0]+"'"
3338        if '2/m' in SytSym and '2' in m:
3339            MagSytSym += '/'
3340        if '-3/m' in SytSym:
3341            MagSytSym = '-'+MagSytSym
3342       
3343    MagSytSym += axis
3344# some exceptions & special rules         
3345    if MagSytSym == "4'/m'm'm'": MagSytSym = "4/m'm'm'"
3346    return MagSytSym
3347   
3348#    if len(GenSym) == 3:
3349#        if SGSpin[1] < 0:
3350#            if 'mm2' in SytSym:
3351#                MagSytSym = "m'm'2"+'('+SplitSytSym[1]
3352#            else:   #bad rule for I41/a
3353#                MagSytSym = SplitSytSym[0]+"'"
3354#                if len(SplitSytSym) > 1:
3355#                    MagSytSym += '('+SplitSytSym[1]
3356#        else:
3357#            MagSytSym = SytSym
3358#        if len(SplitSytSym) >1:
3359#            if "-4'"+'('+SplitSytSym[1] in dupDir:
3360#                MagSytSym = MagSytSym.replace('-4',"-4'")
3361#            if "-6'"+'('+SplitSytSym[1] in dupDir:
3362#                MagSytSym = MagSytSym.replace('-6',"-6'")
3363#        return MagSytSym
3364#           
3365    return SytSym
3366
3367def UpdateSytSym(Phase):
3368    ''' Update site symmetry/site multiplicity after space group/BNS lattice change
3369    '''
3370    generalData = Phase['General']
3371    SGData = generalData['SGData']
3372    Atoms = Phase['Atoms']
3373    cx,ct,cs,cia = generalData['AtomPtrs']
3374    for atom in Atoms:
3375        XYZ = atom[cx:cx+3]
3376        sytsym,Mult = SytSym(XYZ,SGData)[:2]
3377        sytSym,Mul,Nop,dupDir = SytSym(XYZ,SGData)
3378        atom[cs] = sytsym
3379        if generalData['Type'] == 'magnetic':
3380            magSytSym = MagSytSym(sytSym,dupDir,SGData)
3381            atom[cs] = magSytSym
3382        atom[cs+1] = Mult
3383    return
3384   
3385def ElemPosition(SGData):
3386    ''' Under development.
3387    Object here is to return a list of symmetry element types and locations suitable
3388    for say drawing them.
3389    So far I have the element type... getting all possible locations without lookup may be impossible!
3390    '''
3391    Inv = SGData['SGInv']
3392    eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
3393    # get operators & expand if centrosymmetric
3394    Ops = SGData['SGOps']
3395    opM = np.array([op[0].T for op in Ops])
3396    opT = np.array([op[1] for op in Ops])
3397    if Inv:
3398        opM = np.concatenate((opM,-opM))
3399        opT = np.concatenate((opT,-opT))
3400    opMT = list(zip(opM,opT))
3401    for M,T in opMT[1:]:        #skip I
3402        Dt = int(nl.det(M))
3403        Tr = int(np.trace(M))
3404        Dt = -(Dt-1)//2
3405        Es = eleSym[Tr][Dt]
3406        if Dt:              #rotation-inversion
3407            I = np.eye(3)
3408            if Tr == 1:     #mirrors/glides
3409                if np.any(T):       #glide
3410                    M2 = np.inner(M,M)
3411                    MT = np.inner(M,T)+T
3412                    print ('glide',Es,MT)
3413                    print (M2)
3414                else:               #mirror
3415                    print ('mirror',Es,T)
3416                    print (I-M)
3417                X = [-1,-1,-1]
3418            elif Tr == -3:  # pure inversion
3419                X = np.inner(nl.inv(I-M),T)
3420                print ('inversion',Es,X)
3421            else:           #other rotation-inversion
3422                M2 = np.inner(M,M)
3423                MT = np.inner(M,T)+T
3424                print ('rot-inv',Es,MT)
3425                print (M2)
3426                X = [-1,-1,-1]
3427        else:               #rotations
3428            print ('rotation',Es)
3429            X = [-1,-1,-1]
3430        #SymElements.append([Es,X])
3431       
3432    return #SymElements
3433   
3434def ApplyStringOps(A,SGData,X,Uij=[]):
3435    'Needs a doc string'
3436    SGOps = SGData['SGOps']
3437    SGCen = SGData['SGCen']
3438    Ax = A.split('+')
3439    Ax[0] = int(Ax[0])
3440    iC = 1
3441    if Ax[0] < 0:
3442        iC = -1
3443    Ax[0] = abs(Ax[0])
3444    nA = Ax[0]%100-1
3445    cA = Ax[0]//100
3446    Cen = SGCen[cA]
3447    M,T = SGOps[nA]
3448    if len(Ax)>1:
3449        cellA = Ax[1].split(',')
3450        cellA = np.array([int(a) for a in cellA])
3451    else:
3452        cellA = np.zeros(3)
3453    newX = Cen+iC*(np.inner(M,X).T+T)+cellA
3454    if len(Uij):
3455        U = Uij2U(Uij)
3456        U = np.inner(M,np.inner(U,M).T)
3457        newUij = U2Uij(U)
3458        return [newX,newUij]
3459    else:
3460        return newX
3461       
3462def ApplyStringOpsMom(A,SGData,Mom):
3463    'Needs a doc string'
3464    SGOps = SGData['SGOps']
3465    Ax = A.split('+')
3466    Ax[0] = int(Ax[0])
3467    iAx = abs(Ax[0])
3468    nA = iAx%100-1
3469    nC = len(SGOps)*(iAx//100)
3470    NA = nA
3471    if Ax[0] < 0:
3472        NA += len(SGOps)
3473    M,T = SGOps[nA]
3474    if SGData['SGGray']:        #no nonzero moments for gray groups!
3475        newMom = [0.,0.,0.]
3476    else:
3477        newMom = np.inner(Mom,M).T*nl.det(M)*SGData['SpnFlp'][NA+nC]
3478#        print(len(SGOps),Ax[0],iAx,nC,nA,NA,SGData['SpnFlp'][NA],Mom,newMom)
3479    return newMom
3480       
3481def StringOpsProd(A,B,SGData):
3482    """
3483    Find A*B where A & B are in strings '-' + '100*c+n' + '+ijk'
3484    where '-' indicates inversion, c(>0) is the cell centering operator,
3485    n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3486    Should return resultant string - C. SGData - dictionary using entries:
3487
3488       *  'SGCen': cell centering vectors [0,0,0] at least
3489       *  'SGOps': symmetry operations as [M,T] so that M*x+T = x'
3490
3491    """
3492    SGOps = SGData['SGOps']
3493    SGCen = SGData['SGCen']
3494    #1st split out the cell translation part & work on the operator parts
3495    Ax = A.split('+'); Bx = B.split('+')
3496    Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
3497    iC = 0
3498    if Ax[0]*Bx[0] < 0:
3499        iC = 1
3500    Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
3501    nA = Ax[0]%100-1;  nB = Bx[0]%100-1
3502    cA = Ax[0]//100;  cB = Bx[0]//100
3503    Cen = (SGCen[cA]+SGCen[cB])%1.0
3504    cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
3505    Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
3506    Mc = np.inner(Ma,Mb.T)
3507#    print Ma,Mb,Mc
3508    Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
3509#    print Ta,Tb,Tc
3510#    print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
3511    nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
3512    #now the cell translation part
3513    if len(Ax)>1:
3514        cellA = Ax[1].split(',')
3515        cellA = [int(a) for a in cellA]
3516    else:
3517        cellA = [0,0,0]
3518    if len(Bx)>1:
3519        cellB = Bx[1].split(',')
3520        cellB = [int(b) for b in cellB]
3521    else:
3522        cellB = [0,0,0]
3523    cellC = np.add(cellA,cellB)
3524    C = str(((nC+1)+(100*cC))*(1-2*iC))+'+'+ \
3525        str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
3526    return C
3527           
3528def U2Uij(U):
3529    #returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
3530    return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
3531   
3532def Uij2U(Uij):
3533    #returns the thermal motion tensor U from Uij as numpy array
3534    return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
3535
3536def StandardizeSpcName(spcgroup):
3537    '''Accept a spacegroup name where spaces may have not been used
3538    in the names according to the GSAS convention (spaces between symmetry
3539    for each axis) and return the space group name as used in GSAS
3540    '''
3541    rspc = spcgroup.replace(' ','').upper()
3542    # deal with rhombohedral and hexagonal setting designations
3543    rhomb = ''
3544    if rspc[-1:] == 'R':
3545        rspc = rspc[:-1]
3546        rhomb = ' R'
3547    gray = ''
3548    if "1'" in rspc:
3549        gray = " 1'"
3550        rspc = rspc.replace("1'",'')
3551    rspc = rspc.replace("'",'')
3552    if rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
3553        rspc = rspc[:-1]
3554    if rspc[1:3] in ['M3','N3','A3','D3']:      #fix cubic old style
3555        rspc.replace('3','-3')
3556    bns = -1
3557    try:
3558        bns = rspc.index('_')
3559        rspc = rspc.replace(rspc[bns:bns+2],'')
3560    except ValueError:
3561        pass
3562    # look for a match in the spacegroup lists
3563    for i in spglist.values():
3564        for spc in i:
3565            if rspc == spc.replace(' ','').upper():
3566                return spc+gray+rhomb
3567    # how about the post-2002 orthorhombic names?
3568    if rspc in sgequiv_2002_orthorhombic:
3569        return sgequiv_2002_orthorhombic[rspc]+gray
3570    else:
3571    # not found
3572        return ''
3573   
3574def SpaceGroupNumber(spcgroup):
3575    SGNo = -1
3576    SpcGp = StandardizeSpcName(spcgroup)
3577    if not SpcGp:
3578        return SGNo
3579    try:
3580        SGNo = spgbyNum.index(SpcGp)
3581    except ValueError:
3582        pass
3583    return SGNo
3584
3585
3586spgbyNum = []
3587'''Space groups indexed by number'''
3588spgbyNum = [None,
3589        'P 1','P -1',                                                   #1-2
3590        'P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
3591        'C 2/m','P 2/c','P 21/c','C 2/c',                               #3-15
3592        'P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21',
3593        'C 2 2 21','C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
3594        'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
3595        'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
3596        'C m m 2','C m c 21','C c c 2',
3597        'A m m 2','A b m 2','A m a 2','A b a 2',
3598        'F m m 2','F d d 2','I m m 2','I b a 2','I m a 2',
3599        'P m m m','P n n n','P c c m','P b a n',
3600        'P m m a','P n n a','P m n a','P c c a','P b a m','P c c n',
3601        'P b c m','P n n m','P m m n','P b c n','P b c a','P n m a',
3602        'C m c m','C m c a','C m m m','C c c m','C m m a','C c c a',
3603        'F m m m', 'F d d d',
3604        'I m m m','I b a m','I b c a','I m m a',                        #16-74
3605        'P 4','P 41','P 42','P 43',
3606        'I 4','I 41',
3607        'P -4','I -4','P 4/m','P 42/m','P 4/n','P 42/n',
3608        'I 4/m','I 41/a',
3609        'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3610        'P 42 21 2','P 43 2 2','P 43 21 2',
3611        'I 4 2 2','I 41 2 2',
3612        'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c','P 4 n c',
3613        'P 42 m c','P 42 b c',
3614        'I 4 m m','I 4 c m','I 41 m d','I 41 c d',
3615        'P -4 2 m','P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2',
3616        'P -4 c 2','P -4 b 2','P -4 n 2',
3617        'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d',
3618        'P 4/m m m','P 4/m c c','P 4/n b m','P 4/n n c','P 4/m b m',
3619        'P 4/m n c','P 4/n m m','P 4/n c c','P 42/m m c','P 42/m c m',
3620        'P 42/n b c','P 42/n n m','P 42/m b c','P 42/m n m','P 42/n m c',
3621        'P 42/n c m',
3622        'I 4/m m m','I 4/m c m','I 41/a m d','I 41/a c d',
3623        'P 3','P 31','P 32','R 3','P -3','R -3',
3624        'P 3 1 2','P 3 2 1','P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1',
3625        'R 3 2',
3626        'P 3 m 1','P 3 1 m','P 3 c 1','P 3 1 c',
3627        'R 3 m','R 3 c',
3628        'P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1',
3629        'R -3 m','R -3 c',                                               #75-167
3630        'P 6','P 61',
3631        'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3632        'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3633        'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3634        'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',   #168-194
3635        'P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m 3','P n 3',
3636        'F m -3','F d -3','I m -3',
3637        'P a -3','I a -3','P 4 3 2','P 42 3 2','F 4 3 2','F 41 3 2',
3638        'I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2','P -4 3 m',
3639        'F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
3640        'P m -3 m','P n -3 n','P m -3 n','P n -3 m',
3641        'F m -3 m','F m -3 c','F d -3 m','F d -3 c',
3642        'I m -3 m','I a -3 d',]                                       #195-230
3643altSettingOrtho = {}
3644''' A dictionary of alternate settings for orthorhombic unit cells
3645'''
3646altSettingOrtho = {
3647        'P 2 2 2' :{'abc':'P 2 2 2','cab':'P 2 2 2','bca':'P 2 2 2','acb':'P 2 2 2','bac':'P 2 2 2','cba':'P 2 2 2'},
3648        'P 2 2 21' :{'ab