source: trunk/GSASIIspc.py @ 3621

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

Modify Transform tool to allow setting spins & displaying operators
Transform operation matches mag space group selection from Bilbao
select nonstandard orthorhombic works (except P 2221 & P 2'2'21) - get same set of mag atom site symmetries

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