source: trunk/GSASIIspc.py @ 3575

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

add table of alternate space group settings for orthorhombics
remove setting 1--> setting 2 correction for k_SUBGROUPSMAG results; already setting 2.

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