source: trunk/GSASIIspc.py @ 3577

Last change on this file since 3577 was 3577, checked in by vondreele, 3 years ago

delete MAXMAGN - not using it any more
fix (maybe) mag constraints; works for orthogonal axes
fix bug in ReadPowerInstprm?
spelling of Levenberg!
fixes to mag structure from Bilbao stuff
set dmin to be min 1A for UnitCell? show reflections
fix Latt2text to show '-' & Trans2text
fix ext test in GetHistogramPhaseData? for mag structures

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