source: trunk/GSASIIspc.py @ 3591

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

work on FindNonstandard? for mag space groups
enhance UseMagAtomDialog? to show mag site symmetry

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