source: trunk/GSASIIspc.py @ 3703

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

refactoring of GetSSfxuinel

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