source: trunk/GSASIIspc.py @ 3625

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

fix 2 bugs in mag. space group stuff
correct Mag Structure tutorial name

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 170.2 KB
Line 
1# -*- coding: utf-8 -*-
2"""
3*GSASIIspc: Space group module*
4-------------------------------
5
6Space group interpretation routines. Note that space group information is
7stored in a :ref:`Space Group (SGData)<SGData_table>` object.
8
9"""
10########### SVN repository information ###################
11# $Date: 2018-09-27 17:30:16 +0000 (Thu, 27 Sep 2018) $
12# $Author: vondreele $
13# $Revision: 3625 $
14# $URL: trunk/GSASIIspc.py $
15# $Id: GSASIIspc.py 3625 2018-09-27 17:30:16Z 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: 3625 $")
27
28npsind = lambda x: np.sin(x*np.pi/180.)
29npcosd = lambda x: np.cos(x*np.pi/180.)
30nxs = np.newaxis
31twopi = 2.0*np.pi
32DEBUG = False
33   
34################################################################################
35#### Space group codes
36################################################################################
37
38def SpcGroup(SGSymbol):
39    """
40    Determines cell and symmetry information from a short H-M space group name
41
42    :param SGSymbol: space group symbol (string) with spaces between axial fields
43    :returns: (SGError,SGData)
44   
45       * SGError = 0 for no errors; >0 for errors (see SGErrors below for details)
46       * SGData - is a dict (see :ref:`Space Group object<SGData_table>`) with entries:
47       
48             * 'SpGrp': space group symbol, slightly cleaned up
49             * 'SGFixed': True if space group data can not be changed, e.g. from magnetic cif; otherwise False
50             * 'SGGray': True if 1' in symbol - gray group for mag. incommensurate phases
51             * 'SGLaue':  one of '-1', '2/m', 'mmm', '4/m', '4/mmm', '3R',
52               '3mR', '3', '3m1', '31m', '6/m', '6/mmm', 'm3', 'm3m'
53             * 'SGInv': boolean; True if centrosymmetric, False if not
54             * 'SGLatt': one of 'P', 'A', 'B', 'C', 'I', 'F', 'R'
55             * 'SGUniq': one of 'a', 'b', 'c' if monoclinic, '' otherwise
56             * 'SGCen': cell centering vectors [0,0,0] at least
57             * 'SGOps': symmetry operations as [M,T] so that M*x+T = x'
58             * 'SGSys': one of 'triclinic', 'monoclinic', 'orthorhombic',
59               'tetragonal', 'rhombohedral', 'trigonal', 'hexagonal', 'cubic'
60             * 'SGPolax': one of ' ', 'x', 'y', 'x y', 'z', 'x z', 'y z',
61               'xyz', '111' for arbitrary axes
62             * 'SGPtGrp': one of 32 point group symbols (with some permutations), which
63                is filled by SGPtGroup, is external (KE) part of supersymmetry point group
64             * 'SSGKl': default internal (Kl) part of supersymmetry point group; modified
65                in supersymmetry stuff depending on chosen modulation vector for Mono & Ortho
66             * 'BNSlattsym': BNS lattice symbol & cenering op - used for magnetic structures
67
68    """
69    LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
70    LattSym = ('P','A','B','C','I','F','R')
71    UniqSym = ('','','a','b','c','',)
72    SysSym = ('triclinic','monoclinic','orthorhombic','tetragonal','rhombohedral','trigonal','hexagonal','cubic')
73    SGData = {}
74    if len(SGSymbol.split()) < 2:
75        return SGErrors(0),SGData
76    if ':R' in SGSymbol:
77        SGSymbol = SGSymbol.replace(':',' ')    #get rid of ':' in R space group symbols from some cif files
78    SGData['SGGray'] = False
79    if "1'" in SGSymbol:        #set for incommensurate magnetic
80        SGData['SGGray'] = True
81        SGSymbol = SGSymbol.replace("1'",'')
82    SGSymbol = SGSymbol.split(':')[0]   #remove :1/2 setting symbol from some cif files
83    if '-2' in SGSymbol:    #replace bad but legal symbols with correct equivalents
84        SGSymbol = SGSymbol.replace('-2','m')
85    if SGSymbol.split()[1] =='3/m':
86        SGSymbol = SGSymbol.replace('3/m','-6')
87    import pyspg
88    SGInfo = pyspg.sgforpy(SGSymbol)
89    SGData['SpGrp'] = SGSymbol.strip().lower().capitalize()
90    SGData['SGLaue'] = LaueSym[SGInfo[0]-1]
91    SGData['SGInv'] = bool(SGInfo[1])
92    SGData['SGLatt'] = LattSym[SGInfo[2]-1]
93    SGData['SGUniq'] = UniqSym[SGInfo[3]+1]
94    SGData['SGFixed'] = False
95    SGData['SGOps'] = []
96    SGData['SGGen'] = []
97    for i in range(SGInfo[5]):
98        Mat = np.array(SGInfo[6][i])
99        Trns = np.array(SGInfo[7][i])
100        SGData['SGOps'].append([Mat,Trns])
101        if 'array' in str(type(SGInfo[8])):        #patch for old fortran bin?
102            SGData['SGGen'].append(int(SGInfo[8][i]))
103    SGData['BNSlattsym'] = [LattSym[SGInfo[2]-1],[0,0,0]]
104    lattSpin = []
105    if SGData['SGLatt'] == 'P':
106        SGData['SGCen'] = np.array(([0,0,0],))
107    elif SGData['SGLatt'] == 'A':
108        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5]))
109        lattSpin += [1,]
110    elif SGData['SGLatt'] == 'B':
111        SGData['SGCen'] = np.array(([0,0,0],[.5,0,.5]))
112        lattSpin += [1,]
113    elif SGData['SGLatt'] == 'C':
114        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,0,]))
115        lattSpin += [1,]
116    elif SGData['SGLatt'] == 'I':
117        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,.5]))
118        lattSpin += [1,]
119    elif SGData['SGLatt'] == 'F':
120        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5],[.5,0,.5],[.5,.5,0,]))
121        lattSpin += [1,1,1,1]
122    elif SGData['SGLatt'] == 'R':
123        SGData['SGCen'] = np.array(([0,0,0],[2./3,1./3,1./3],[1./3,2./3,2./3]))
124
125    if SGData['SGInv']:
126        if SGData['SGLaue'] in ['-1','2/m','mmm']:
127            Ibar = 7
128        elif SGData['SGLaue'] in ['4/m','4/mmm']:
129            Ibar = 1
130        elif SGData['SGLaue'] in ['3R','3mR','3','3m1','31m','6/m','6/mmm']:
131            Ibar = 15 #8+4+2+1
132        else:
133            Ibar = 4
134        Ibarx = Ibar&14
135    else:
136        Ibarx = 8
137        if SGData['SGLaue'] in ['-1','2/m','mmm','m3','m3m']:
138            Ibarx = 0
139    moregen = []
140    for i,gen in enumerate(SGData['SGGen']):
141        if SGData['SGLaue'] in ['m3','m3m']:
142            if gen in [1,2,4]:
143                SGData['SGGen'][i] = 4
144            elif gen < 7:
145                SGData['SGGen'][i] = 0
146        elif SGData['SGLaue'] in ['4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm']:
147            if gen == 2:
148                SGData['SGGen'][i] = 4
149            elif gen in [3,5]:
150                SGData['SGGen'][i] = 3
151            elif gen == 6:
152                if SGData['SGLaue'] in ['4/m','4/mmm']:
153                    SGData['SGGen'][i] = 128
154                else:
155                    SGData['SGGen'][i] = 16
156            elif not SGData['SGInv'] and gen == 12:
157                SGData['SGGen'][i] = 8
158            elif (not SGData['SGInv']) and (SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']) and (gen == 1):
159                SGData['SGGen'][i] = 24
160        gen = SGData['SGGen'][i]
161        if gen == 99:
162            gen = 8
163            if SGData['SGLaue'] in ['3m1','31m','6/m','6/mmm']:
164                gen = 3
165            elif SGData['SGLaue'] == 'm3m':
166                gen = 12
167            SGData['SGGen'][i] = gen
168        elif gen == 98:
169            gen = 8
170            if SGData['SGLaue'] in ['3m1','31m','6/m','6/mmm']:
171                gen = 4
172            SGData['SGGen'][i] = gen
173        elif not SGData['SGInv'] and gen in [23,] and SGData['SGLaue'] in ['m3','m3m']:
174            SGData['SGGen'][i] = 24
175        elif gen >= 16 and gen != 128:
176            if not SGData['SGInv']:
177                gen = 31
178            else:
179                gen ^= Ibarx
180            SGData['SGGen'][i] = gen
181        if SGData['SGInv']:
182            if gen < 128:
183                moregen.append(SGData['SGGen'][i]^Ibar)
184            else:
185                moregen.append(1)
186    SGData['SGGen'] += moregen
187    if SGData['SGLaue'] in '-1':
188        SGData['SGSys'] = SysSym[0]
189    elif SGData['SGLaue'] in '2/m':
190        SGData['SGSys'] = SysSym[1]
191    elif SGData['SGLaue'] in 'mmm':
192        SGData['SGSys'] = SysSym[2]
193    elif SGData['SGLaue'] in ['4/m','4/mmm']:
194        SGData['SGSys'] = SysSym[3]
195    elif SGData['SGLaue'] in ['3R','3mR']:
196        SGData['SGSys'] = SysSym[4]
197    elif SGData['SGLaue'] in ['3','3m1','31m']:
198        SGData['SGSys'] = SysSym[5]
199    elif SGData['SGLaue'] in ['6/m','6/mmm']:
200        SGData['SGSys'] = SysSym[6]
201    elif SGData['SGLaue'] in ['m3','m3m']:
202        SGData['SGSys'] = SysSym[7]
203    SGData['SGPolax'] = SGpolar(SGData)
204    SGData['SGPtGrp'],SGData['SSGKl'] = SGPtGroup(SGData)
205
206    if SGData['SGLatt'] == 'R':
207        if SGData['SGPtGrp'] in ['3',]:
208            SGData['SGSpin'] = 3*[1,]
209        elif SGData['SGPtGrp'] in  ['-3','32','3m']:
210            SGData['SGSpin'] = 4*[1,]
211        elif SGData['SGPtGrp'] in  ['-3m',]:
212            SGData['SGSpin'] = 5*[1,]
213       
214    else:
215        if SGData['SGPtGrp'] in ['1','3','23',]:
216            SGData['SGSpin'] = lattSpin+[1,]
217        elif SGData['SGPtGrp'] in ['-1','2','m','4','-4','-3','312','321','3m1','31m','6','-6','432','-43m']:
218            SGData['SGSpin'] = lattSpin+[1,1,]
219        elif SGData['SGPtGrp'] in ['2/m','4/m','422','4mm','-42m','-4m2','-3m1','-31m',
220            '6/m','622','6mm','-6m2','-62m','m3','m3m']:
221            SGData['SGSpin'] = lattSpin+[1,1,1,]
222        else: #'222'-'mmm','4/mmm','6/mmm'
223            SGData['SGSpin'] = lattSpin+[1,1,1,1,]
224    return SGInfo[-1],SGData
225
226def SGErrors(IErr):
227    '''
228    Interprets the error message code from SpcGroup. Used in SpaceGroup.
229   
230    :param IErr: see SGError in :func:`SpcGroup`
231    :returns:
232        ErrString - a string with the error message or "Unknown error"
233    '''
234
235    ErrString = [' ',
236        'Less than 2 operator fields were found',
237        'Illegal Lattice type, not P, A, B, C, I, F or R',
238        'Rhombohedral lattice requires a 3-axis',
239        'Minus sign does not preceed 1, 2, 3, 4 or 6',
240        'Either a 5-axis anywhere or a 3-axis in field not allowed',
241        ' ',
242        'I for COMPUTED GO TO out of range.',
243        'An a-glide mirror normal to A not allowed',
244        'A b-glide mirror normal to B not allowed',
245        'A c-glide mirror normal to C not allowed',
246        'D-glide in a primitive lattice not allowed',
247        'A 4-axis not allowed in the 2nd operator field',
248        'A 6-axis not allowed in the 2nd operator field',
249        'More than 24 matrices needed to define group',
250        ' ',
251        'Improper construction of a rotation operator',
252        'Mirror following a / not allowed',
253        'A translation conflict between operators',
254        'The 2bar operator is not allowed',
255        '3 fields are legal only in R & m3 cubic groups',
256        'Syntax error. Expected I -4 3 d at this point',
257        ' ',
258        'A or B centered tetragonal not allowed',
259        ' ','unknown error in sgroup',' ',' ',' ',
260        'Illegal character in the space group symbol',
261        ]
262    try:
263        return ErrString[IErr]
264    except:
265        return "Unknown error"
266
267def SGpolar(SGData):
268    '''
269    Determine identity of polar axes if any
270    '''
271    POL = ('','x','y','x y','z','x z','y z','xyz','111')
272    NP = [1,2,4]
273    NPZ = [0,1]
274    for M,T in SGData['SGOps']:
275        for i in range(3):
276            if M[i][i] <= 0.: NP[i] = 0
277        if M[0][2] > 0: NPZ[0] = 8
278        if M[1][2] > 0: NPZ[1] = 0
279    NPol = (NP[0]+NP[1]+NP[2]+NPZ[0]*NPZ[1])*(1-int(SGData['SGInv']))
280    return POL[NPol]
281   
282def SGPtGroup(SGData):
283    '''
284    Determine point group of the space group - done after space group symbol has
285    been evaluated by SpcGroup. Only short symbols are allowed
286   
287    :param SGData: from :func SpcGroup
288    :returns: SSGPtGrp & SSGKl (only defaults for Mono & Ortho)
289    '''
290    Flds = SGData['SpGrp'].split()
291    if len(Flds) < 2:
292        return '',[]
293    if SGData['SGLaue'] == '-1':    #triclinic
294        if '-' in Flds[1]:
295            return '-1',[-1,]
296        else:
297            return '1',[1,]
298    elif SGData['SGLaue'] == '2/m': #monoclinic - default for 2D modulation vector
299        if '/' in SGData['SpGrp']:
300            return '2/m',[-1,1]
301        elif '2' in SGData['SpGrp']:
302            return '2',[-1,]
303        else:
304            return 'm',[1,]
305    elif SGData['SGLaue'] == 'mmm': #orthorhombic
306        if SGData['SpGrp'].count('2') == 3:
307            return '222',[-1,-1,-1]
308        elif SGData['SpGrp'].count('2') == 1:
309            if SGData['SGPolax'] == 'x':
310                return '2mm',[-1,1,1]
311            elif SGData['SGPolax'] == 'y':
312                return 'm2m',[1,-1,1]
313            elif SGData['SGPolax'] == 'z':
314                return 'mm2',[1,1,-1]
315        else:
316            return 'mmm',[1,1,1]
317    elif SGData['SGLaue'] == '4/m': #tetragonal
318        if '/' in SGData['SpGrp']:
319            return '4/m',[1,-1]
320        elif '-' in Flds[1]:
321            return '-4',[-1,]
322        else:
323            return '4',[1,]
324    elif SGData['SGLaue'] == '4/mmm':
325        if '/' in SGData['SpGrp']:
326            return '4/mmm',[1,-1,1,1]
327        elif '-' in Flds[1]:
328            if '2' in Flds[2]:
329                return '-42m',[-1,-1,1]
330            else:
331                return '-4m2',[-1,1,-1]             
332        elif '2' in Flds[2:]:
333            return '422',[1,-1,-1]
334        else:
335            return '4mm',[1,1,1]
336    elif SGData['SGLaue'] in ['3','3R']:  #trigonal/rhombohedral
337        if '-' in Flds[1]:
338            return '-3',[-1,]
339        else:
340            return '3',[1,]
341    elif SGData['SGLaue'] == '3mR' or 'R' in Flds[0]:
342        if '2' in Flds[2]:
343            return '32',[1,-1]
344        elif '-' in Flds[1]:
345            return '-3m',[-1,1]
346        else:
347            return '3m',[1,1]
348    elif SGData['SGLaue'] == '3m1':
349        if '2' in Flds[2]:
350            return '321',[1,-1,1]
351        elif '-' in Flds[1]:
352            return '-3m1',[-1,1,1]
353        else:
354            return '3m1',[1,1,1]
355    elif SGData['SGLaue'] == '31m':
356        if '2' in Flds[3]:
357            return '312',[1,1,-1]
358        elif '-' in Flds[1]:
359            return '-31m',[-1,1,1]
360        else:
361            return '31m',[1,1,1]
362    elif SGData['SGLaue'] == '6/m': #hexagonal
363        if '/' in SGData['SpGrp']:
364            return '6/m',[1,-1]
365        elif '-' in SGData['SpGrp']:
366            return '-6',[-1,]
367        else:
368            return '6',[1,]
369    elif SGData['SGLaue'] == '6/mmm':
370        if '/' in SGData['SpGrp']:
371            return '6/mmm',[1,-1,1,1]
372        elif '-' in Flds[1]:
373            if '2' in Flds[2]:
374                return '-62m',[-1,-1,1]
375            else:
376                return '-6m2',[-1,1,-1]                 
377        elif '2' in Flds[2:]:
378            return '622',[1,-1,-1]
379        else:
380            return '6mm',[1,1,1]   
381    elif SGData['SGLaue'] == 'm3':      #cubic - no (3+1) supersymmetry
382        if '2' in Flds[1]:
383            return '23',[]
384        else: 
385            return 'm3',[]
386    elif SGData['SGLaue'] == 'm3m':
387        if '4' in Flds[1]:
388            if '-' in Flds[1]:
389                return '-43m',[]
390            else:
391                return '432',[]
392        else:
393            return 'm3m',[]
394   
395def SGPrint(SGData,AddInv=False):
396    '''
397    Print the output of SpcGroup in a nicely formatted way. Used in SpaceGroup
398
399    :param SGData: from :func:`SpcGroup`
400    :returns:
401        SGText - list of strings with the space group details
402        SGTable - list of strings for each of the operations
403    '''
404    if SGData.get('SGFixed',False):       #inverses included in ops for cif fixed
405        Mult = len(SGData['SGCen'])*len(SGData['SGOps'])
406    else:
407        Mult = len(SGData['SGCen'])*len(SGData['SGOps'])*(int(SGData['SGInv'])+1)
408    SGText = []
409    SGText.append(' Space Group: '+SGData['SpGrp'])
410    if SGData.get('SGGray',False): SGText[-1] += " 1'"
411    CentStr = 'centrosymmetric'
412    if not SGData['SGInv']:
413        CentStr = 'non'+CentStr
414    if SGData['SGLatt'] in 'ABCIFR':
415        SGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
416    else:
417        SGText.append(' The lattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower()) 
418    SGText.append(' The Laue symmetry is '+SGData['SGLaue'])
419    if 'SGPtGrp' in SGData:         #patch
420        SGText.append(' The lattice point group is '+SGData['SGPtGrp'])
421    SGText.append(' Multiplicity of a general site is '+str(Mult))
422    if SGData['SGUniq'] in ['a','b','c']:
423        SGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
424    if SGData['SGInv']:
425        SGText.append(' The inversion center is located at 0,0,0')
426    if SGData['SGPolax']:
427        SGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
428    SGText.append(' ')
429    if len(SGData['SGCen']) == 1:
430        SGText.append(' The equivalent positions are:\n')
431    else:   
432        SGText.append(' The equivalent positions are:\n')
433        SGText.append(' ('+Latt2text(SGData['SGCen'])+')+\n')
434    SGTable = []
435    for i,Opr in enumerate(SGData['SGOps']):
436        SGTable.append('(%2d) %s'%(i+1,MT2text(Opr)))
437    if AddInv and SGData['SGInv']:
438        for i,Opr in enumerate(SGData['SGOps']):
439            IOpr = [-Opr[0],-Opr[1]]
440            SGTable.append('(%2d) %s'%(i+1,MT2text(IOpr)))       
441#    if SGData.get('SGGray',False) and not SGData.get('SGFixed',False):
442#        SGTable.append("     for 1'")
443#        for i,Opr in enumerate(SGData['SGOps']):
444#            SGTable.append('(%2d) %s'%(i+1,MT2text(Opr)))
445#        if AddInv and SGData['SGInv']:
446#            for i,Opr in enumerate(SGData['SGOps']):
447#                IOpr = [-Opr[0],-Opr[1]]
448#                SGTable.append('(%2d) %s'%(i+1,MT2text(IOpr)))       
449    return SGText,SGTable
450
451def AllOps(SGData):
452    '''
453    Returns a list of all operators for a space group, including those for
454    centering and a center of symmetry
455   
456    :param SGData: from :func:`SpcGroup`
457    :returns: (SGTextList,offsetList,symOpList,G2oprList) where
458
459      * SGTextList: a list of strings with formatted and normalized
460        symmetry operators.
461      * offsetList: a tuple of (dx,dy,dz) offsets that relate the GSAS-II
462        symmetry operation to the operator in SGTextList and symOpList.
463        these dx (etc.) values are added to the GSAS-II generated
464        positions to provide the positions that are generated
465        by the normalized symmetry operators.       
466      * symOpList: a list of tuples with the normalized symmetry
467        operations as (M,T) values
468        (see ``SGOps`` in the :ref:`Space Group object<SGData_table>`)
469      * G2oprList: The GSAS-II operations for each symmetry operation as
470        a tuple with (center,mult,opnum,opcode), where center is (0,0,0), (0.5,0,0),
471        (0.5,0.5,0.5),...; where mult is 1 or -1 for the center of symmetry
472        where opnum is the number for the symmetry operation, in ``SGOps``
473        (starting with 0) and opcode is mult*(100*icen+j+1).
474    '''
475    SGTextList = []
476    offsetList = []
477    symOpList = []
478    G2oprList = []
479    G2opcodes = []
480    onebar = (1,)
481    if SGData['SGInv']:
482        onebar += (-1,)
483    for icen,cen in enumerate(SGData['SGCen']):
484        for mult in onebar:
485            for j,(M,T) in enumerate(SGData['SGOps']):
486                offset = [0,0,0]
487                Tprime = (mult*T)+cen
488                for i in range(3):
489                    while Tprime[i] < 0:
490                        Tprime[i] += 1
491                        offset[i] += 1
492                    while Tprime[i] >= 1:
493                        Tprime[i] += -1
494                        offset[i] += -1
495                Opr = [mult*M,Tprime]
496                OPtxt = MT2text(Opr)
497                SGTextList.append(OPtxt.replace(' ',''))
498                offsetList.append(tuple(offset))
499                symOpList.append((mult*M,Tprime))
500                G2oprList.append((cen,mult,j))
501                G2opcodes.append(mult*(100*icen+j+1))
502    return SGTextList,offsetList,symOpList,G2oprList,G2opcodes
503
504def TextOps(text,table,reverse=False):
505    ''' Makes formatted operator list
506        :param text,table: arrays of text made by SGPrint
507        :param reverse: True for x+1/2 form; False for 1/2+x form
508        :returns: OpText: full list of symmetry operators; one operation per line
509        generally printed to console for use via cut/paste in other programs, but
510        could be used for direct input
511    '''
512    OpText = []
513    Inv = True
514    if 'noncentro' in text[1]:
515        Inv = False
516    Cent = [[0,0,0],]
517    if '0,0,0' in text[-1]:
518        Cent = np.array(eval(text[-1].split('+')[0].replace(';','),(')))
519    OpsM = []
520    OpsT = []
521    for item in table:
522        M,T = Text2MT(item.split(')')[1].replace(' ',''),CIF=True)
523        OpsM.append(M)
524        OpsT.append(T)
525    OpsM = np.array(OpsM)
526    OpsT = np.array(OpsT)
527    if Inv:
528        OpsM = np.concatenate((OpsM,-OpsM))
529        OpsT = np.concatenate((OpsT,-OpsT%1.))
530    for cent in Cent:
531        for iop,opM in enumerate(list(OpsM)):
532            txt = MT2text([opM,(OpsT[iop]+cent)%1.],reverse)
533            OpText.append(txt.replace(' ','').lower())
534    return OpText
535
536def TextGen(SGData,reverse=False):
537    GenSym,GenFlg,BNSsym = GetGenSym(SGData)
538    SGData['GenSym'] = GenSym
539    SGData['GenFlg'] = GenFlg
540    text,table = SGPrint(SGData)
541    GenText = []
542    OprNames = GetOprNames(SGData)
543    OpText = TextOps(text,table,reverse)
544    for name in SGData['GenSym']:
545        gid = OprNames.index(name.replace(' ',''))
546        GenText.append(OpText[gid])
547    if len(SGData['SGCen']) > 1:
548        GenText.append(OpText[-1])
549    return GenText
550
551def GetOprNames(SGData):
552    OprNames = [GetOprPtrName(str(irtx)) for irtx in PackRot(SGData['SGOps'])]
553    if SGData['SGInv']:
554        OprNames += [GetOprPtrName(str(-irtx)) for irtx in PackRot(SGData['SGOps'])]
555    return OprNames
556   
557def MT2text(Opr,reverse=False):
558    "From space group matrix/translation operator returns text version"
559    XYZ = ('-Z','-Y','-X','X-Y','ERR','Y-X','X','Y','Z')
560    TRA = ('   ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
561    Fld = ''
562    M,T = Opr
563    for j in range(3):
564        IJ = int(round(2*M[j][0]+3*M[j][1]+4*M[j][2]+4))%12
565        IK = int(round(T[j]*12))%12
566        if IK:
567            if IJ < 3:
568                if reverse:
569                    Fld += (XYZ[IJ]+'+'+TRA[IK]).rjust(5)
570                else:
571                    Fld += (TRA[IK]+XYZ[IJ]).rjust(5)
572            else:
573                if reverse:
574                    Fld += (XYZ[IJ]+'+'+TRA[IK]).rjust(5)
575                else:
576                    Fld += (TRA[IK]+'+'+XYZ[IJ]).rjust(5)
577        else:
578            Fld += XYZ[IJ].rjust(5)
579        if j != 2: Fld += ', '
580    return Fld
581   
582def Latt2text(Cen):
583    "From lattice centering vectors returns ';' delimited cell centering vectors"
584    lattTxt = ''
585    fracList = ['1/2','1/3','2/3','1/4','3/4','1/5','2/5','3/5','4/5','1/6','5/6',
586        '1/7','2/7','3/7','4/7','5/7','6/7','1/8','3/8','5/8','7/8','1/9','2/9','4/9','5/9','7/9','8/9']
587    mulList = [2,3,3,4,4,5,5,5,5,6,6,7,7,7,7,7,7,8,8,8,8,9,9,9,9,9,9]
588    prodList = [1.,1.,2.,1.,3.,1.,2.,3.,4.,1.,5.,1.,2.,3.,4.,5.,6.,1.,3.,5.,7.,1.,2.,4.,5.,7.,8.]
589    nCen = len(Cen)
590    for i,cen in enumerate(Cen):
591        txt = ''
592        for icen in cen:
593            if icen == 1:
594                txt += '1,'
595                continue
596            if not icen:
597                txt += '0,'
598                continue
599            if icen < 0:
600                txt += '-'
601                icen *= -1
602            for mul,prod,frac in zip(mulList,prodList,fracList):
603                if abs(icen*mul-prod) < 1.e-5:
604                    txt += frac+','
605                    break
606        lattTxt += txt[:-1]+'; '
607        if i and not i%8 and i < nCen-1:    #not for the last cen!
608            lattTxt += '\n     '
609    return lattTxt[:-2]
610
611def SpaceGroup(SGSymbol):
612    '''
613    Print the output of SpcGroup in a nicely formatted way.
614
615    :param SGSymbol: space group symbol (string) with spaces between axial fields
616    :returns: nothing
617    '''
618    E,A = SpcGroup(SGSymbol)
619    if E > 0:
620        print (SGErrors(E))
621        return
622    for l in SGPrint(A):
623        print (l)
624################################################################################
625#### Magnetic space group stuff
626################################################################################
627       
628def SetMagnetic(SGData):
629    GenSym,GenFlg,BNSsym = GetGenSym(SGData)
630    SGData['GenSym'] = GenSym
631    SGData['GenFlg'] = GenFlg
632    OprNames,SpnFlp = GenMagOps(SGData)
633    SGData['SpnFlp'] = SpnFlp
634    SGData['MagSpGrp'] = MagSGSym(SGData)
635
636def GetGenSym(SGData):
637    '''
638    Get the space group generator symbols
639    :param SGData: from :func:`SpcGroup`
640    LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
641    LattSym = ('P','A','B','C','I','F','R')
642   
643    '''
644    OprNames = [GetOprPtrName(str(irtx)) for irtx in PackRot(SGData['SGOps'])]
645    if SGData['SGInv']:
646        OprNames += [GetOprPtrName(str(-irtx)) for irtx in PackRot(SGData['SGOps'])]
647    Nsyms = len(SGData['SGOps'])
648    if SGData['SGInv'] and not SGData['SGFixed']: Nsyms *= 2
649    UsymOp = ['1',]
650    OprFlg = [0,] 
651    if Nsyms == 2:                    #Centric triclinic or acentric monoclinic
652        UsymOp.append(OprNames[1])
653        OprFlg.append(SGData['SGGen'][1])
654    elif Nsyms == 4:                    #Point symmetry 2/m, 222, 22m, or 4
655        if '4z' in OprNames[1]:          #Point symmetry 4 or -4
656            UsymOp.append(OprNames[1])
657            OprFlg.append(SGData['SGGen'][1])
658        elif not SGData['SGInv']:       #Acentric Orthorhombic
659            if 'm' in OprNames[1:4]:    #22m, 2m2 or m22
660                if '2' in OprNames[1]:      #Acentric orthorhombic, 2mm
661                    UsymOp.append(OprNames[2])
662                    OprFlg.append(SGData['SGGen'][2])
663                    UsymOp.append(OprNames[3])
664                    OprFlg.append(SGData['SGGen'][3])
665                elif '2' in OprNames[2]:    #Acentric orthorhombic, m2m
666                    UsymOp.append(OprNames[1])
667                    OprFlg.append(SGData['SGGen'][1])
668                    UsymOp.append(OprNames[3])
669                    OprFlg.append(SGData['SGGen'][3])
670                else:                       #Acentric orthorhombic, mm2
671                    UsymOp.append(OprNames[1])
672                    OprFlg.append(SGData['SGGen'][1])
673                    UsymOp.append(OprNames[2])
674                    OprFlg.append(SGData['SGGen'][2])
675            else:                           #Acentric orthorhombic, 222
676                SGData['SGGen'][1:] = [4,2,1]
677                UsymOp.append(OprNames[1])
678                OprFlg.append(SGData['SGGen'][1])
679                UsymOp.append(OprNames[2])
680                OprFlg.append(SGData['SGGen'][2])
681                UsymOp.append(OprNames[3])
682                OprFlg.append(SGData['SGGen'][3])
683        else:                               #Centric Monoclinic
684            UsymOp.append(OprNames[1])
685            OprFlg.append(SGData['SGGen'][1])
686            UsymOp.append(OprNames[3])
687            OprFlg.append(SGData['SGGen'][3])
688    elif Nsyms == 6:                    #Point symmetry 32, 3m or 6
689            if '6' in OprNames[1]:      #Hexagonal 6/m Laue symmetry
690                UsymOp.append(OprNames[1])
691                OprFlg.append(SGData['SGGen'][1])
692            else:                       #Trigonal
693                UsymOp.append(OprNames[4])
694                OprFlg.append(SGData['SGGen'][3])
695                if '2110' in OprNames[1]: UsymOp[-1] = ' 2100 '
696    elif Nsyms == 8:                    #Point symmetry mmm, 4/m, or 422, etc
697        if '4' in OprNames[1]:           #Tetragonal
698            if SGData['SGInv']:         #4/m
699                UsymOp.append(OprNames[1])
700                OprFlg.append(SGData['SGGen'][1])
701                UsymOp.append(OprNames[6])
702                OprFlg.append(SGData['SGGen'][6])
703            else:
704                if 'x' in OprNames[4]:      #4mm type group
705                    UsymOp.append(OprNames[4])
706                    OprFlg.append(6)
707                    UsymOp.append(OprNames[7])
708                    OprFlg.append(8)
709                else:                       #-42m, -4m2, and 422 type groups
710                    UsymOp.append(OprNames[5])
711                    OprFlg.append(8)
712                    UsymOp.append(OprNames[6])
713                    OprFlg.append(19)
714        else:                               #Orthorhombic, mmm
715            UsymOp.append(OprNames[1])
716            OprFlg.append(SGData['SGGen'][1])
717            UsymOp.append(OprNames[2])
718            OprFlg.append(SGData['SGGen'][2])
719            UsymOp.append(OprNames[7])
720            OprFlg.append(SGData['SGGen'][7])
721    elif Nsyms == 12 and '3' in OprNames[1] and SGData['SGSys'] != 'cubic':        #Trigonal
722        UsymOp.append(OprNames[3])
723        OprFlg.append(SGData['SGGen'][3])
724        UsymOp.append(OprNames[9])
725        OprFlg.append(SGData['SGGen'][9])
726    elif Nsyms == 12 and '6' in OprNames[1]:        #Hexagonal
727        if 'mz' in OprNames[9]:                     #6/m
728            UsymOp.append(OprNames[1])
729            OprFlg.append(SGData['SGGen'][1])
730            UsymOp.append(OprNames[6])
731            OprFlg.append(SGData['SGGen'][6])
732        else:                                       #6mm, -62m, -6m2 or 622
733            UsymOp.append(OprNames[6])
734            OprFlg.append(18)
735            if 'm' in UsymOp[-1]: OprFlg[-1] = 20
736            UsymOp.append(OprNames[7])
737            OprFlg.append(24)
738    elif Nsyms in [16,24]:
739        if '3' in OprNames[1]:
740            UsymOp.append('')
741            OprFlg.append(SGData['SGGen'][3])
742            for i in range(Nsyms):
743                if 'mx' in OprNames[i]:
744                    UsymOp[-1] = OprNames[i]
745                elif 'm11' in OprNames[i]:
746                    UsymOp[-1] = OprNames[i]
747                elif '211' in OprNames[i]:
748                    UsymOp[-1] = OprNames[i]
749                    OprFlg[-1] = 24
750        else:                                     #4/mmm or 6/mmm
751            UsymOp.append('  mz  ')
752            OprFlg.append(1)
753            if '4' in OprNames[1]:                  #4/mmm
754                UsymOp.append('  mx  ')
755                OprFlg.append(20)
756                UsymOp.append(' m110 ')
757                OprFlg.append(24)
758            else:                                   #6/mmm
759                UsymOp.append(' m110 ')
760                OprFlg.append(4)
761                UsymOp.append(' m+-0 ')
762                OprFlg.append(8)
763    else:                                           #System is cubic
764        if Nsyms == 48:
765            UsymOp.append('  mx  ')
766            OprFlg.append(4)
767            UsymOp.append(' m110 ')
768            OprFlg.append(24)
769           
770    if 'P' in SGData['SGLatt']:
771        if SGData['SGSys'] == 'triclinic':
772            BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5]}           
773        elif SGData['SGSys'] == 'monoclinic':
774            BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5]}
775            if SGData['SGUniq'] == 'a':
776                BNSsym.update({'P_B':[.5,0,.5],'P_C':[.5,.5,0]})
777            elif SGData['SGUniq'] == 'b':
778                BNSsym.update({'P_A':[.5,.5,0],'P_C':[0,.5,.5]})
779            elif SGData['SGUniq'] == 'c':
780                BNSsym.update({'P_A':[0,.5,.5],'P_B':[.5,0,.5]})
781        elif SGData['SGSys'] == 'orthorhombic':
782            BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5],
783                'P_A':[0,.5,.5],'P_B':[.5,0,.5],'P_C':[.5,.5,0],'P_I':[.5,.5,.5]}
784        elif SGData['SGSys'] == 'tetragonal':
785            BNSsym = {'P_c':[0,0,.5],'P_C':[.5,.5,0],'P_I':[.5,.5,.5]}           
786        elif SGData['SGSys'] in ['trigonal','hexagonal']:
787            BNSsym = {'P_c':[0,0,.5]}           
788        elif SGData['SGSys'] == 'cubic':
789            BNSsym = {'P_I':[.5,.5,.5]}           
790           
791    elif 'A' in SGData['SGLatt']:
792        if SGData['SGSys'] == 'monoclinic':
793            BNSsym = {}
794            if SGData['SGUniq'] == 'b':
795                BNSsym.update({'A_a':[.5,0,0],'A_c':[0,0,.5]})
796            elif SGData['SGUniq'] == 'c':
797                BNSsym.update({'A_a':[.5,0,0],'A_b':[0,.5,0]})
798        elif SGData['SGSys'] == 'orthorhombic':
799            BNSsym = {'A_a':[.5,0,0],'A_b':[0,.5,0],'A_c':[0,0,.5],
800               'A_B':[.5,0,.5],'A_C':[.5,.5,0]}   
801        elif SGData['SGSys'] == 'triclinic':
802            BNSsym = {'A_a':[.5,0,0],'A_b':[0,.5,0],'A_c':[0,0,.5]}   
803           
804    elif 'B' in SGData['SGLatt']:
805        if SGData['SGSys'] == 'monoclinic':
806            BNSsym = {}
807            if SGData['SGUniq'] == 'a':
808                BNSsym.update({'B_b':[0,.5,0],'B_c':[0,0,.5]})
809            elif SGData['SGUniq'] == 'c':
810                BNSsym.update({'B_a':[.5,0,0],'B_b':[0,.5,0]})
811        elif SGData['SGSys'] == 'orthorhombic':
812            BNSsym = {'B_a':[.5,0,0],'B_b':[0,.5,0],'B_c':[0,0,.5],
813                'B_A':[0,.5,.5],'B_C':[.5,.5,0]}     
814        elif SGData['SGSys'] == 'triclinic':
815            BNSsym = {'B_a':[.5,0,0],'B_b':[0,.5,0],'B_c':[0,0,.5]}     
816           
817    elif 'C' in SGData['SGLatt']:
818        if SGData['SGSys'] == 'monoclinic':
819            BNSsym = {}
820            if SGData['SGUniq'] == 'a':
821                BNSsym.update({'C_b':[0,.5,.0],'C_c':[0,0,.5]})
822            elif SGData['SGUniq'] == 'b':
823                BNSsym.update({'C_a':[.5,0,0],'C_c':[0,0,.5]})
824        elif SGData['SGSys'] == 'orthorhombic':
825            BNSsym = {'C_a':[.5,0,0],'C_b':[0,.5,0],'C_c':[0,0,.5],
826                'C_A':[0,.5,.5],'C_B':[.5,0,.5]}     
827        elif SGData['SGSys'] == 'triclinic':
828            BNSsym = {'C_a':[.5,0,0],'C_b':[0,.5,0],'C_c':[0,0,.5]}     
829           
830    elif 'I' in SGData['SGLatt']:
831        if SGData['SGSys'] in ['monoclinic','orthorhombic','triclinic']:
832            BNSsym = {'I_a':[.5,0,0],'I_b':[0,.5,0],'I_c':[0,0,.5]}
833        elif SGData['SGSys'] == 'tetragonal':
834            BNSsym = {'I_c':[0,0,.5]}
835        elif SGData['SGSys'] == 'cubic':
836            BNSsym = {} 
837           
838    elif 'F' in SGData['SGLatt']:
839        if SGData['SGSys'] in ['monoclinic','orthorhombic','cubic','triclinic']:
840            BNSsym = {'F_S':[.5,.5,.5]}
841           
842    elif 'R' in SGData['SGLatt']:
843        BNSsym = {'R_I':[0,0,.5]}
844    return UsymOp,OprFlg,BNSsym
845
846def ApplyBNSlatt(SGData,BNSlatt):
847    Tmat = np.eye(3)
848    BNS = BNSlatt[0]
849    A = np.array(BNSlatt[1])
850    Laue = SGData['SGLaue']
851    SGCen = SGData['SGCen']
852    if '_a' in BNS:
853        Tmat[0,0] = 2.0
854    elif '_b' in BNS:
855        Tmat[1,1] = 2.0
856    elif '_c' in BNS:
857        Tmat[2,2] = 2.0
858    elif '_A' in BNS:
859        Tmat[0,0] = 2.0
860    elif '_B' in BNS:
861        Tmat[1,1] = 2.0
862    elif '_C' in BNS:
863        Tmat[2,2] = 2.0
864    elif '_I' in BNS:
865        Tmat *= 2.0
866        if 'R' in Laue:
867            SGData['SGSpin'][-1] = -1
868        else:
869            SGData['SGSpin'].append(-1)
870    elif '_S' in BNS:
871        SGData['SGSpin'][-1] = -1
872        SGData['SGSpin'] += [-1,-1,-1,]
873        Tmat *= 2.0
874    else:
875        return Tmat
876    SGData['SGSpin'].append(-1)
877    C = SGCen+A
878    SGData['SGCen'] = np.vstack((SGCen,C))%1.
879    return Tmat
880       
881def CheckSpin(isym,SGData):
882    ''' Check for exceptions in spin rules
883    '''
884    if SGData['SGPtGrp'] in ['222','mm2','2mm','m2m']:      #only 2/3 can be red; not 1/3 or 3/3
885        if SGData['SGSpin'][1]*SGData['SGSpin'][2]*SGData['SGSpin'][3] < 0:
886            SGData['SGSpin'][(isym+1)%3+1] *= -1
887        if SGData['SpGrp'][0] == 'F' and isym > 2:
888            SGData['SGSpin'][(isym+1)%3+3] == 1
889    elif SGData['SGPtGrp'] == 'mmm':
890        if SGData['SpGrp'][0] == 'F' and isym > 2:
891            SGData['SGSpin'][(isym+1)%3+3] == 1
892
893def MagSGSym(SGData):       #needs to use SGPtGrp not SGLaue!
894    SGLaue = SGData['SGLaue']
895    if '1' not in SGData['GenSym']:        #patch for old gpx files
896        SGData['GenSym'] = ['1',]+SGData['GenSym']
897        SGData['SGSpin'] = [1,]+list(SGData['SGSpin'])
898    if len(SGData['SGSpin'])<len(SGData['GenSym']):
899        SGData['SGSpin'] = [1,]+list(SGData['SGSpin'])      #end patch
900    GenSym = SGData['GenSym'][1:]       #skip identity
901    SpnFlp = SGData['SGSpin']
902#    print('SpnFlp',SpnFlp)
903    SGPtGrp = SGData['SGPtGrp']
904    if len(SpnFlp) == 1:
905        SGData['MagPtGp'] = SGPtGrp
906        return SGData['SpGrp']
907    magSym = SGData['SpGrp'].split()
908    if SGLaue in ['-1',]:
909        SGData['MagPtGp'] = SGPtGrp
910        if SpnFlp[1] == -1:
911            magSym[1] += "'"
912            SGData['MagPtGp'] += "'"
913    elif SGLaue in ['2/m','4/m','6/m']: #all ok
914        Uniq = {'a':1,'b':2,'c':3,'':1}
915        id = [0,1]
916        if len(magSym) > 2:
917            id = [0,Uniq[SGData['SGUniq']]]
918        sym = magSym[id[1]].split('/')
919        Ptsym = SGLaue.split('/')
920        if len(GenSym) == 3:
921            for i in [0,1,2]:
922                if SpnFlp[i+1] < 0:
923                    sym[i] += "'"
924                    Ptsym[i] += "'"
925        else:
926            for i in range(len(GenSym)):
927                if SpnFlp[i+1] < 0:                     
928                    sym[i] += "'"
929                    Ptsym[i] += "'"
930        SGData['MagPtGp'] = '/'.join(Ptsym)
931        magSym[id[1]] = '/'.join(sym)
932    elif SGPtGrp in ['mmm','mm2','m2m','2mm','222']:
933        SGData['MagPtGp'] = ''
934        for i in [0,1,2]:
935            SGData['MagPtGp'] += SGPtGrp[i]
936            if SpnFlp[i+1] < 0:
937                magSym[i+1] += "'"
938                SGData['MagPtGp'] += "'"
939    elif SGLaue == '6/mmm': #ok
940        magPtGp = list(SGPtGrp)
941        if len(GenSym) == 2:
942            for i in [0,1]:
943                if SpnFlp[i+1] < 0:
944                    magSym[i+2] += "'"
945                    magPtGp[i+1] += "'"
946            if SpnFlp[1]*SpnFlp[2] < 0:
947                magSym[1] += "'"
948                magPtGp[0] += "'"
949        else:
950            sym = magSym[1].split('/')
951            Ptsym = ['6','m']
952            magPtGp = ['','m','m']
953            for i in [0,1,2]:
954                if SpnFlp[i+1] < 0:
955                    if i:
956                        magSym[i+1] += "'"
957                        magPtGp[i] += "'"
958                    else:
959                        sym[1] += "'"
960                        Ptsym[0] += "'"
961            if SpnFlp[2]*SpnFlp[3] < 0:
962                sym[0] += "'"                   
963                Ptsym[0] += "'"                   
964            magSym[1] = '/'.join(sym)
965            magPtGp[0] = '/'.join(Ptsym)
966        SGData['MagPtGp'] = ''.join(magPtGp)
967    elif SGLaue == '4/mmm':
968        magPtGp = list(SGPtGrp)
969        if len(GenSym) == 2:
970            for i in [0,1]:
971                if SpnFlp[i+1] < 0:
972                    magSym[i+2] += "'"
973                    magPtGp[i+1] += "'"
974            if SpnFlp[1]*SpnFlp[2] < 0:
975                magSym[1] += "'"
976                magPtGp[0] += "'"
977        else:
978            if '/' in magSym[1]:    #P 4/m m m, etc.
979                sym = magSym[1].split('/')
980                Ptsym = ['4','m']
981                magPtGp = ['','m','m']
982                for i in [0,1,2]:
983                    if SpnFlp[i+1] < 0:
984                        if i:
985                            magSym[i+1] += "'"
986                            magPtGp[i] += "'"
987                        else:
988                            sym[1] += "'"
989                            Ptsym[1] += "'"
990                if SpnFlp[2]*SpnFlp[3] < 0:
991                    sym[0] += "'"                   
992                    Ptsym[0] += "'"                   
993                magSym[1] = '/'.join(sym)
994                magPtGp[0] = '/'.join(Ptsym)
995            else:
996                for i in [0,1]:
997                    if SpnFlp[i+1] < 0:
998                        magSym[i+2] += "'"
999                if SpnFlp[1]*SpnFlp[2] < 0:
1000                    magSym[1] += "'"
1001        SGData['MagPtGp'] = ''.join(magPtGp)
1002    elif SGLaue in ['3','3m1','31m']:   #ok
1003        Ptsym = list(SGLaue)
1004        if len(GenSym) == 1:    #all ok
1005            id = 2
1006            if (len(magSym) == 4) and (magSym[2] == '1'):
1007                id = 3
1008            if '3' in GenSym[0]:
1009                id = 1
1010            magSym[id].strip("'")
1011            if SpnFlp[1] < 0:
1012                magSym[id] += "'"
1013                Ptsym[id-1] += "'"
1014        elif len(GenSym) == 2:
1015            if 'R' in GenSym[1]:
1016                magSym[-1].strip("'")
1017                if SpnFlp[1] < 0:
1018                    magSym[-1] += "'"
1019                    Ptsym[-1] += "'"
1020            else:
1021                i,j = [1,2]
1022                if magSym[2] == '1':
1023                    i,j = [1,3]
1024                magSym[i].strip("'")
1025                Ptsym[i-1].strip("'")
1026                magSym[j].strip("'")
1027                Ptsym[j-1].strip("'")
1028                if SpnFlp[1:3] == [1,-1]:
1029                    magSym[i] += "'"
1030                    Ptsym[i-1] += "'"
1031                elif SpnFlp[1:3] == [-1,-1]:
1032                    magSym[j] += "'"
1033                    Ptsym[j-1] += "'"
1034                elif SpnFlp[1:3] == [-1,1]:
1035                    magSym[i] += "'"
1036                    Ptsym[i-1] += "'"
1037                    magSym[j] += "'"
1038                    Ptsym[j-1] += "'"
1039        elif len(GenSym):
1040            if 'c' not in magSym[2]:
1041                i,j = [1,2]
1042                magSym[i].strip("'")
1043                Ptsym[i-1].strip("'")
1044                magSym[j].strip("'")
1045                Ptsym[j-1].strip("'")
1046                if SpnFlp[1:3] == [1,-1]:
1047                    magSym[i] += "'"
1048                    Ptsym[i-1] += "'"
1049                elif SpnFlp[1:3] == [-1,-1]:
1050                    magSym[j] += "'"
1051                    Ptsym[j-1] += "'"
1052                elif SpnFlp[2] == [-1,1]:
1053                    magSym[i] += "'"
1054                    Ptsym[i-1] += "'"
1055                    magSym[j] += "'"
1056                    Ptsym[j-1] += "'"
1057        SGData['MagPtGp'] = ''.join(Ptsym)
1058    elif SGData['SGPtGrp'] == '23' and len(magSym):
1059        SGData['MagPtGp'] = '23'
1060    elif SGData['SGPtGrp'] == 'm3':
1061        SGData['MagPtGp'] = "m3"
1062        if SpnFlp[1] < 0:
1063            magSym[1] += "'"
1064            magSym[2] += "'"
1065            SGData['MagPtGp'] = "m'3'"
1066        if SpnFlp[2] < 0:
1067            if not 'm' in magSym[1]:    #only Ia3
1068                magSym[1].strip("'")
1069                SGData['MagPtGp'] = "m3'"
1070    elif SGData['SGPtGrp'] in ['432','-43m']:
1071        Ptsym = SGData['SGPtGrp'].split('3')
1072        if SpnFlp[1] < 0:
1073            magSym[1] += "'"
1074            Ptsym[0] += "'"
1075            magSym[3] += "'"
1076            Ptsym[1] += "'"
1077        SGData['MagPtGp'] = '3'.join(Ptsym)
1078    elif SGData['SGPtGrp'] == 'm3m':
1079        Ptsym = ['m','3','m']
1080        if SpnFlp[1:3] == [-1,1]:
1081            magSym[1] += "'"
1082            Ptsym[0] += "'"
1083            magSym[2] += "'"
1084            Ptsym[1] += "'"
1085        elif SpnFlp[1:3] == [1,-1]:
1086            magSym[3] += "'"
1087            Ptsym[2] += "'"
1088        elif SpnFlp[1:3] == [-1,-1]:
1089            magSym[1] += "'"
1090            Ptsym[0] += "'"
1091            magSym[2] += "'"
1092            Ptsym[1] += "'"
1093            magSym[3] += "'"
1094            Ptsym[2] += "'"
1095        SGData['MagPtGp'] = ''.join(Ptsym)
1096#    print SpnFlp
1097    magSym[0] = SGData.get('BNSlattsym',[SGData['SGLatt'],[0,0,0]])[0]
1098    return ' '.join(magSym)
1099
1100def fixMono(SpGrp):
1101    'fixes b-unique monoclinics in e.g. P 1 2/1c 1 --> P 21/c '
1102    Flds = SpGrp.split()
1103    if len(Flds) == 4:
1104        if Flds[2] != '1':
1105            return '%s %s'%(Flds[0],Flds[2])
1106        else:
1107            return None
1108    else:
1109        return SpGrp
1110
1111def Trans2Text(Trans):
1112    "from transformation matrix to text"
1113    cells = ['a','b','c']
1114    Text = ''
1115    for row in Trans:
1116        Fld = ''
1117        for i in [0,1,2]:
1118            if row[i]:
1119                if Fld and row[i] > 0.:
1120                    Fld += '+'
1121                Fld += '%3.1f'%(row[i])+cells[i]
1122        Text += Fld
1123        Text += ','
1124        Text = Text.replace('1.0','').replace('.0','').replace('0.5','1/2')
1125    return Text[:-1]
1126
1127def getlattSym(Trans):
1128    Fives = {'ababc':'abc','bcbca':'cba','acacb':'acb','cabab':'cab','abcab':'acb'}
1129    transText = Trans2Text(Trans)
1130    lattSym = ''
1131    for fld in transText.split(','):
1132        if 'a' in fld: lattSym += 'a'
1133        if 'b' in fld: lattSym += 'b'
1134        if 'c' in fld: lattSym += 'c'
1135    if len(lattSym) != 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        dF = None
2560        dFTP = None
2561        if siteSym == '1':
2562            CSI = [[1,0],[2,0]],2*[[1.,0.],]
2563        elif siteSym == '-1':
2564            CSI = [[1,0],[0,0]],2*[[1.,0.],]
2565        else:
2566            delt2 = np.eye(2)*0.001
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        dX,dXTP = None,None
2618        if siteSym == '1':
2619            CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2620        elif siteSym == '-1':
2621            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.],]
2622        else:
2623            delt4 = np.ones(4)*0.001
2624            delt5 = np.ones(5)*0.001
2625            delt6 = np.eye(6)*0.001
2626            if 'Fourier' in waveType:
2627                dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2628                  #3x6x12 modulated position array (X,Spos,tau)& force positive
2629                CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2630            elif waveType == 'Sawtooth':
2631                dX = posSawtooth(tau,delt4[0],delt4[1:])
2632                CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0]]),
2633                    np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2634            elif waveType in ['ZigZag','Block']:
2635                if waveType == 'ZigZag':
2636                    dX = posZigZag(tau,delt5[:2],delt5[2:])
2637                else:
2638                    dX = posBlock(tau,delt5[:2],delt5[2:])
2639                CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2640                    np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2641            XSC = np.ones(6,dtype='i')
2642            dXTP = []
2643            for i in SdIndx:
2644                sop = Sop[i]
2645                ssop = SSop[i]
2646                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2647                xsc = np.ones(6,dtype='i')
2648                if 'Fourier' in waveType:
2649                    dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,nxs,nxs]
2650                elif waveType == 'Sawtooth':
2651                    dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,nxs,nxs]
2652                elif waveType == 'ZigZag':
2653                    dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2654                elif waveType == 'Block':
2655                    dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2656                dXT = np.inner(sop[0],dXT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2657                dXT = np.swapaxes(dXT,1,2)      # back to array(3x6x49)
2658                dXT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2659                dXTP.append(dXT)
2660                if waveType == 'Fourier':
2661                    for i in range(3):
2662                        if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2663                            xsc[i] = 0
2664                        if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2665                            xsc[i+3] = 0
2666                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2667                        xsc[3:6] = 0
2668                        CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2669                            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2670                        if dT:
2671                            if '(x)' in siteSym:
2672                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2673                                if 'm' in siteSym and len(SdIndx) == 1:
2674                                    CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2675                            elif '(y)' in siteSym:
2676                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2677                                if 'm' in siteSym and len(SdIndx) == 1:
2678                                    CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2679                            elif '(z)' in siteSym:
2680                                CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2681                                if 'm' in siteSym and len(SdIndx) == 1:
2682                                    CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2683                        else:
2684                            CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2685                    if '4/mmm' in laue:
2686                        if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2687                            if '(xy)' in siteSym:
2688                                CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2689                                if dT:
2690                                    CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2691                                else:
2692                                    CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2693                        if '(xy)' in siteSym or '(+-0)' in siteSym:
2694                            mul = 1
2695                            if '(+-0)' in siteSym:
2696                                mul = -1
2697                            if np.allclose(dX[0,0,:],dXT[1,0,:]):
2698                                CSI[0][3:5] = [[11,0,0],[11,0,0]]
2699                                CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2700                                xsc[3:5] = 0
2701                            if np.allclose(dX[0,3,:],dXT[0,4,:]):
2702                                CSI[0][:2] = [[12,0,0],[12,0,0]]
2703                                CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2704                                xsc[:2] = 0
2705                XSC &= xsc
2706                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2707            if waveType == 'Fourier':
2708                n = -1
2709                if debug: print (XSC)
2710                for i,X in enumerate(XSC):
2711                    if X:
2712                        n += 1
2713                        CSI[0][i][0] = n+1
2714                        CSI[1][i][0] = 1.0
2715           
2716        return list(CSI),dX,dXTP
2717       
2718    def DoUij():
2719        dU,dUTP = None,None
2720        if siteSym == '1':
2721            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2722                [7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2723        elif siteSym == '-1':
2724            CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],   \
2725                6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2726        else:
2727            tau = np.linspace(0,1,49,True)
2728            delt12 = np.eye(12)*0.0001
2729            dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
2730            CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2731            USC = np.ones(12,dtype='i')
2732            dUTP = []
2733            dtau = 0.
2734            for i in SdIndx:
2735                sop = Sop[i]
2736                ssop = SSop[i]
2737                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2738                usc = np.ones(12,dtype='i')
2739                dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
2740                dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3)    #convert dUT to 12x49x3x3
2741                dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2742                dUT = np.array(U2Uij(dUijT))    #convert to 6x12x49
2743                dUT = dUT[:,:,np.argsort(tauT)]
2744                dUT[:,:6,:] *=(ssdet*sdet)
2745                dUTP.append(dUT)
2746                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2747                    if dT:
2748                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2749                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2750                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2751                        [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2752                    else:
2753                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2754                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2755                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2756                        [0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2757                    if 'mm2(x)' in siteSym and dT:
2758                        CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2759                        USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2760                    elif '(xy)' in siteSym and dT:
2761                        CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2762                            [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2763                        CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2764                        USC = [1,1,1,1,1,1,1,1,1,1,1,1]                             
2765                    elif '(x)' in siteSym and dT:
2766                        CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2767                    elif '(y)' in siteSym and dT:
2768                        CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2769                    elif '(z)' in siteSym and dT:
2770                        CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2771                    for i in range(6):
2772                        if not USC[i]:
2773                            CSI[0][i] = [0,0,0]
2774                            CSI[1][i] = [0.,0.,0.]
2775                            CSI[0][i+6] = [0,0,0]
2776                            CSI[1][i+6] = [0.,0.,0.]
2777                else:                       
2778                    for i in range(6):
2779                        if not np.allclose(dU[i,i,:],dUT[i,i,:]):  #sin part
2780                            usc[i] = 0
2781                        if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]):   #cos part
2782                            usc[i+6] = 0
2783                    if np.any(dUT[1,0,:]):
2784                        if '4/m' in siteSym:
2785                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2786                            if ssop[1][3]:
2787                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2788                                usc[9] = 1
2789                            else:
2790                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2791                                usc[9] = 0
2792                        elif '4' in siteSym:
2793                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2794                            CSI[0][:2] = [[11,0,0],[11,0,0]]
2795                            if ssop[1][3]:
2796                                CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2797                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2798                                usc[2] = 0
2799                                usc[8] = 0
2800                                usc[3] = 1
2801                                usc[9] = 1
2802                            else:
2803                                CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2804                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2805                                usc[2] = 1
2806                                usc[8] = 1
2807                                usc[3] = 0               
2808                                usc[9] = 0
2809                        elif 'xy' in siteSym or '+-0' in siteSym:
2810                            if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2811                                CSI[0][4:6] = [[12,0,0],[12,0,0]]
2812                                CSI[0][6:8] = [[11,0,0],[11,0,0]]
2813                                CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2814                                CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2815                                usc[4:6] = 0
2816                                usc[6:8] = 0
2817                           
2818                    if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2819                USC &= usc
2820            if debug: print (USC)
2821            if not np.any(dtau%.5):
2822                n = -1
2823                for i,U in enumerate(USC):
2824                    if U:
2825                        n += 1
2826                        CSI[0][i][0] = n+1
2827                        CSI[1][i][0] = 1.0
2828   
2829        return list(CSI),dU,dUTP
2830   
2831    def DoMag():
2832        CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2833        if siteSym == '1':
2834            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2835        elif siteSym in ['-1','mmm','2/m(x)','2/m(y)','2/m(z)','4/mmm001']:
2836            CSI = 3*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3*[[0.,0.,0.],]+3*[[1.,0.,0.],]
2837        else:
2838            delt6 = np.eye(6)*0.001
2839            dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,nxs,nxs]
2840              #3x6x12 modulated moment array (M,Spos,tau)& force positive
2841            CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2842            MSC = np.ones(6,dtype='i')
2843            dMTP = []
2844            for i in SdIndx:
2845                sop = Sop[i]
2846                ssop = SSop[i]
2847                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2848                msc = np.ones(6,dtype='i')
2849                dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,nxs,nxs]
2850                dMT = np.inner(sop[0],dMT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2851                dMT = np.swapaxes(dMT,1,2)      # back to array(3x6x49)
2852                dMT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2853                dMTP.append(dMT)
2854                for i in range(3):
2855                    if not np.allclose(dM[i,i,:],dMT[i,i,:]):
2856                        msc[i] = 0
2857                    if not np.allclose(dM[i,i+3,:],dMT[i,i+3,:]):
2858                        msc[i+3] = 0
2859                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2860                    msc[3:6] = 0
2861                    CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2862                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2863                    if dT:
2864                        if '(x)' in siteSym:
2865                            CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2866                            if 'm' in siteSym and len(SdIndx) == 1:
2867                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2868                        elif '(y)' in siteSym:
2869                            CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2870                            if 'm' in siteSym and len(SdIndx) == 1:
2871                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2872                        elif '(z)' in siteSym:
2873                            CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2874                            if 'm' in siteSym and len(SdIndx) == 1:
2875                                CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2876                    else:
2877                        CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2878                if '4/mmm' in laue:
2879                    if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2880                        if '(xy)' in siteSym:
2881                            CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2882                            if dT:
2883                                CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2884                            else:
2885                                CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2886                    if '(xy)' in siteSym or '(+-0)' in siteSym:
2887                        mul = 1
2888                        if '(+-0)' in siteSym:
2889                            mul = -1
2890                        if np.allclose(dM[0,0,:],dMT[1,0,:]):
2891                            CSI[0][3:5] = [[11,0,0],[11,0,0]]
2892                            CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2893                            msc[3:5] = 0
2894                        if np.allclose(dM[0,3,:],dMT[0,4,:]):
2895                            CSI[0][:2] = [[12,0,0],[12,0,0]]
2896                            CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2897                            msc[:2] = 0
2898                MSC &= msc
2899                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2900            n = -1
2901            if debug: print (MSC)
2902            for i,M in enumerate(MSC):
2903                if M:
2904                    n += 1
2905                    CSI[0][i][0] = n+1
2906                    CSI[1][i][0] = 1.0
2907
2908        return CSI,None,None
2909       
2910    if debug: print ('super space group: '+SSGData['SSpGrp'])
2911    xyz = np.array(XYZ)%1.
2912    SGOps = copy.deepcopy(SGData['SGOps'])
2913    laue = SGData['SGLaue']
2914    siteSym = SytSym(XYZ,SGData)[0].strip()
2915    if debug: print ('siteSym: '+siteSym)
2916    SSGOps = copy.deepcopy(SSGData['SSGOps'])
2917    #expand ops to include inversions if any
2918    if SGData['SGInv']:
2919        for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
2920            SGOps.append([-op[0],-op[1]%1.])
2921            SSGOps.append([-sop[0],-sop[1]%1.])
2922    #build set of sym ops around special position       
2923    SSop = []
2924    Sop = []
2925    Sdtau = []
2926    for iop,Op in enumerate(SGOps):         
2927        nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
2928        if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
2929            SSop.append(SSGOps[iop])
2930            Sop.append(SGOps[iop])
2931            ssopinv = nl.inv(SSGOps[iop][0])
2932            mst = ssopinv[3][:3]
2933            epsinv = ssopinv[3][3]
2934            Sdtau.append(np.sum(mst*(XYZ-SGOps[iop][1])-epsinv*SSGOps[iop][1][3]))
2935    SdIndx = np.argsort(np.array(Sdtau))     # just to do in sensible order
2936    if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
2937    #setup displacement arrays
2938    tau = np.linspace(-1,1,49,True)
2939    #make modulation arrays - one parameter at a time
2940    if Stype == 'Sfrac':
2941        CSI,dF,dFTP = DoFrac()
2942    elif Stype == 'Spos':
2943        CSI,dF,dFTP = DoXYZ()
2944        CSI[0] = orderParms(CSI[0])
2945    elif Stype == 'Sadp':
2946        CSI,dF,dFTP = DoUij()
2947        CSI[0] = orderParms(CSI[0]) 
2948    elif Stype == 'Smag':
2949        CSI,dF,dFTP = DoMag()
2950       
2951    if debug:
2952        return CSI,dF,dFTP
2953    else:
2954        return CSI
2955   
2956def MustrainNames(SGData):
2957    'Needs a doc string'
2958    laue = SGData['SGLaue']
2959    uniq = SGData['SGUniq']
2960    if laue in ['m3','m3m']:
2961        return ['S400','S220']
2962    elif laue in ['6/m','6/mmm','3m1']:
2963        return ['S400','S004','S202']
2964    elif laue in ['31m','3']:
2965        return ['S400','S004','S202','S211']
2966    elif laue in ['3R','3mR']:
2967        return ['S400','S220','S310','S211']
2968    elif laue in ['4/m','4/mmm']:
2969        return ['S400','S004','S220','S022']
2970    elif laue in ['mmm']:
2971        return ['S400','S040','S004','S220','S202','S022']
2972    elif laue in ['2/m']:
2973        SHKL = ['S400','S040','S004','S220','S202','S022']
2974        if uniq == 'a':
2975            SHKL += ['S013','S031','S211']
2976        elif uniq == 'b':
2977            SHKL += ['S301','S103','S121']
2978        elif uniq == 'c':
2979            SHKL += ['S130','S310','S112']
2980        return SHKL
2981    else:
2982        SHKL = ['S400','S040','S004','S220','S202','S022']
2983        SHKL += ['S310','S103','S031','S130','S301','S013']
2984        SHKL += ['S211','S121','S112']
2985        return SHKL
2986       
2987def HStrainVals(HSvals,SGData):
2988    laue = SGData['SGLaue']
2989    uniq = SGData['SGUniq']
2990    DIJ = np.zeros(6)
2991    if laue in ['m3','m3m']:
2992        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
2993    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2994        DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
2995    elif laue in ['3R','3mR']:
2996        DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
2997    elif laue in ['4/m','4/mmm']:
2998        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
2999    elif laue in ['mmm']:
3000        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3001    elif laue in ['2/m']:
3002        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3003        if uniq == 'a':
3004            DIJ[5] = HSvals[3]
3005        elif uniq == 'b':
3006            DIJ[4] = HSvals[3]
3007        elif uniq == 'c':
3008            DIJ[3] = HSvals[3]
3009    else:
3010        DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
3011    return DIJ
3012
3013def HStrainNames(SGData):
3014    'Needs a doc string'
3015    laue = SGData['SGLaue']
3016    uniq = SGData['SGUniq']
3017    if laue in ['m3','m3m']:
3018        return ['D11','eA']         #add cubic strain term
3019    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3020        return ['D11','D33']
3021    elif laue in ['3R','3mR']:
3022        return ['D11','D12']
3023    elif laue in ['4/m','4/mmm']:
3024        return ['D11','D33']
3025    elif laue in ['mmm']:
3026        return ['D11','D22','D33']
3027    elif laue in ['2/m']:
3028        Dij = ['D11','D22','D33']
3029        if uniq == 'a':
3030            Dij += ['D23']
3031        elif uniq == 'b':
3032            Dij += ['D13']
3033        elif uniq == 'c':
3034            Dij += ['D12']
3035        return Dij
3036    else:
3037        Dij = ['D11','D22','D33','D12','D13','D23']
3038        return Dij
3039   
3040def MustrainCoeff(HKL,SGData):
3041    'Needs a doc string'
3042    #NB: order of terms is the same as returned by MustrainNames
3043    laue = SGData['SGLaue']
3044    uniq = SGData['SGUniq']
3045    h,k,l = HKL
3046    Strm = []
3047    if laue in ['m3','m3m']:
3048        Strm.append(h**4+k**4+l**4)
3049        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
3050    elif laue in ['6/m','6/mmm','3m1']:
3051        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
3052        Strm.append(l**4)
3053        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
3054    elif laue in ['31m','3']:
3055        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
3056        Strm.append(l**4)
3057        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
3058        Strm.append(4.0*h*k*l*(h+k))
3059    elif laue in ['3R','3mR']:
3060        Strm.append(h**4+k**4+l**4)
3061        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
3062        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
3063        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
3064    elif laue in ['4/m','4/mmm']:
3065        Strm.append(h**4+k**4)
3066        Strm.append(l**4)
3067        Strm.append(3.0*(h*k)**2)
3068        Strm.append(3.0*((h*l)**2+(k*l)**2))
3069    elif laue in ['mmm']:
3070        Strm.append(h**4)
3071        Strm.append(k**4)
3072        Strm.append(l**4)
3073        Strm.append(3.0*(h*k)**2)
3074        Strm.append(3.0*(h*l)**2)
3075        Strm.append(3.0*(k*l)**2)
3076    elif laue in ['2/m']:
3077        Strm.append(h**4)
3078        Strm.append(k**4)
3079        Strm.append(l**4)
3080        Strm.append(3.0*(h*k)**2)
3081        Strm.append(3.0*(h*l)**2)
3082        Strm.append(3.0*(k*l)**2)
3083        if uniq == 'a':
3084            Strm.append(2.0*k*l**3)
3085            Strm.append(2.0*l*k**3)
3086            Strm.append(4.0*k*l*h**2)
3087        elif uniq == 'b':
3088            Strm.append(2.0*l*h**3)
3089            Strm.append(2.0*h*l**3)
3090            Strm.append(4.0*h*l*k**2)
3091        elif uniq == 'c':
3092            Strm.append(2.0*h*k**3)
3093            Strm.append(2.0*k*h**3)
3094            Strm.append(4.0*h*k*l**2)
3095    else:
3096        Strm.append(h**4)
3097        Strm.append(k**4)
3098        Strm.append(l**4)
3099        Strm.append(3.0*(h*k)**2)
3100        Strm.append(3.0*(h*l)**2)
3101        Strm.append(3.0*(k*l)**2)
3102        Strm.append(2.0*k*h**3)
3103        Strm.append(2.0*h*l**3)
3104        Strm.append(2.0*l*k**3)
3105        Strm.append(2.0*h*k**3)
3106        Strm.append(2.0*l*h**3)
3107        Strm.append(2.0*k*l**3)
3108        Strm.append(4.0*k*l*h**2)
3109        Strm.append(4.0*h*l*k**2)
3110        Strm.append(4.0*k*h*l**2)
3111    return Strm
3112
3113def MuShklMean(SGData,Amat,Shkl):
3114   
3115    def genMustrain(xyz,Shkl):
3116        uvw = np.inner(Amat.T,xyz)
3117        Strm = np.array(MustrainCoeff(uvw,SGData))
3118        Sum = np.sum(np.multiply(Shkl,Strm))
3119        Sum = np.where(Sum > 0.01,Sum,0.01)
3120        Sum = np.sqrt(Sum)
3121        return Sum*xyz
3122       
3123    PHI = np.linspace(0.,360.,30,True)
3124    PSI = np.linspace(0.,180.,30,True)
3125    X = np.outer(npcosd(PHI),npsind(PSI))
3126    Y = np.outer(npsind(PHI),npsind(PSI))
3127    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3128    XYZ = np.dstack((X,Y,Z))
3129    XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
3130    return np.sqrt(np.sum(XYZ**2)/900.)
3131   
3132def Muiso2Shkl(muiso,SGData,cell):
3133    "this is to convert isotropic mustrain to generalized Shkls"
3134    import GSASIIlattice as G2lat
3135    A = G2lat.cell2AB(cell)[0]
3136   
3137    def minMus(Shkl,muiso,H,SGData,A):
3138        U = np.inner(A.T,H)
3139        S = np.array(MustrainCoeff(U,SGData))
3140        Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,nxs]),axis=0))
3141        rad = np.sqrt(np.sum((Sum[:,nxs]*H)**2,axis=1))
3142        return (muiso-rad)**2
3143       
3144    laue = SGData['SGLaue']
3145    PHI = np.linspace(0.,360.,60,True)
3146    PSI = np.linspace(0.,180.,60,True)
3147    X = np.outer(npsind(PHI),npsind(PSI))
3148    Y = np.outer(npcosd(PHI),npsind(PSI))
3149    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3150    HKL = np.dstack((X,Y,Z))
3151    if laue in ['m3','m3m']:
3152        S0 = [1000.,1000.]
3153    elif laue in ['6/m','6/mmm','3m1']:
3154        S0 = [1000.,1000.,1000.]
3155    elif laue in ['31m','3']:
3156        S0 = [1000.,1000.,1000.,1000.]
3157    elif laue in ['3R','3mR']:
3158        S0 = [1000.,1000.,1000.,1000.]
3159    elif laue in ['4/m','4/mmm']:
3160        S0 = [1000.,1000.,1000.,1000.]
3161    elif laue in ['mmm']:
3162        S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
3163    elif laue in ['2/m']:
3164        S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
3165    else:
3166        S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000., 
3167            1000.,1000.,0.,0.,0.]
3168    S0 = np.array(S0)
3169    HKL = np.reshape(HKL,(-1,3))
3170    result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
3171    return result[0]
3172       
3173def PackRot(SGOps):
3174    IRT = []
3175    for ops in SGOps:
3176        M = ops[0]
3177        irt = 0
3178        for j in range(2,-1,-1):
3179            for k in range(2,-1,-1):
3180                irt *= 3
3181                irt += M[k][j]
3182        IRT.append(int(irt))
3183    return IRT
3184       
3185def SytSym(XYZ,SGData):
3186    '''
3187    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
3188
3189    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
3190    :param SGData: from SpcGroup
3191    :Returns: a four element tuple:
3192
3193     * The 1st element is a code for the site symmetry (see GetKNsym)
3194     * The 2nd element is the site multiplicity
3195     * Ndup number of overlapping operators
3196     * dupDir Dict - dictionary of overlapping operators
3197
3198    '''
3199    Mult = 1
3200    Isym = 0
3201    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
3202        Isym = 1073741824
3203    Jdup = 0
3204    Ndup = 0
3205    dupDir = {}
3206    inv = SGData['SGInv']+1
3207    icen = SGData['SGCen']
3208    if SGData['SGFixed']:       #already in list of operators
3209        inv = 1
3210    Xeqv = list(GenAtom(XYZ,SGData,True))
3211#    for xeqv in Xeqv:   print(xeqv)
3212    IRT = PackRot(SGData['SGOps'])
3213    L = -1
3214    for ic,cen in enumerate(icen):
3215        for invers in range(int(inv)):
3216            for io,ops in enumerate(SGData['SGOps']):
3217                irtx = (1-2*invers)*IRT[io]
3218                L += 1
3219                if not Xeqv[L][1]:
3220                    Ndup = io
3221                    Jdup += 1
3222                    jx = GetOprPtrNumber(str(irtx))   #[KN table no,op name,KNsym ptr]
3223                    if jx < 39:
3224                        px = GetOprName(str(irtx))
3225                        if Xeqv[L][-1] < 0:
3226                            if '(' in px:
3227                                px = px.split('(')
3228                                px[0] += "'"
3229                                px = '('.join(px)
3230                            else:   
3231                                px += "'"
3232                        dupDir[px] = L
3233                        Isym += 2**(jx-1)
3234    if Isym == 1073741824: Isym = 0
3235    Mult = len(SGData['SGOps'])*len(SGData['SGCen'])*inv//Jdup
3236         
3237    return GetKNsym(str(Isym)),Mult,Ndup,dupDir
3238   
3239def MagSytSym(SytSym,dupDir,SGData):
3240    '''
3241    site sym operations: 1,-1,2,3,-3,4,-4,6,-6,m need to be marked if spin inversion
3242    '''
3243    SGData['GenSym'],SGData['GenFlg'] = GetGenSym(SGData)[:2]
3244#    print('SGPtGrp',SGData['SGPtGrp'],'SytSym',SytSym,'MagSpGrp',SGData['MagSpGrp'])
3245#    print('dupDir',dupDir)
3246    SplitSytSym = SytSym.split('(')
3247    if SGData['SGGray']:
3248        return SytSym+"1'"
3249    if SytSym == '1':       #genersl position
3250        return SytSym
3251    if SplitSytSym[0] == SGData['SGPtGrp']:     #simple cases
3252        try:
3253            MagSytSym = SGData['MagSpGrp'].split()[1]
3254        except IndexError:
3255            MagSytSym = SGData['MagSpGrp'][1:].strip("1'")
3256        if len(SplitSytSym) > 1:
3257            MagSytSym += '('+SplitSytSym[1]
3258        return MagSytSym
3259    if len(dupDir) == 1:
3260        return list(dupDir.keys())[0]
3261   
3262   
3263    if '2/m' in SytSym:         #done I think; last 2wo might be not needed
3264        ops = {'(x)':['2(x)','m(x)'],'(y)':['2(y)','m(y)'],'(z)':['2(z)','m(z)'],
3265               '(100)':['2(100)','m(100)'],'(010)':['2(010)','m(010)'],'(001)':['2(001)','m(001)'],
3266               '(120)':['2(120)','m(120)'],'(210)':['2(210)','m(210)'],'(+-0)':['2(+-0)','m(+-0)'],
3267               '(110)':['2(110)','m(110)']}
3268   
3269    elif '4/mmm' in SytSym:
3270        ops = {'(x)':['4(x)','m(x)','m(y)','m(0+-)'],   #m(0+-) for cubic m3m?
3271               '(y)':['4(y)','m(y)','m(z)','m(+0-)'],   #m(+0-)
3272               '(z)':['4(z)','m(z)','m(x)','m(+-0)']}   #m(+-0)
3273    elif '4mm' in SytSym:
3274        ops = {'(x)':['4(x)','m(y)','m(yz)'],'(y)':['4(y)','m(z)','m(xz)'],'(z)':['4(z)','m(x)','m(110)']}
3275    elif '422' in SytSym:
3276        ops = {'(x)':['4(x)','2(y)','2(yz)'],'(y)':['4(y)','2(z)','2(xz)'],'(z)':['4(z)','2(x)','2(110)']}
3277    elif '-4m2' in SytSym:
3278        ops = {'(x)':['-4(x)','m(x)','2(yz)'],'(y)':['-4(y)','m(y)','2(xz)'],'(z)':['-4(z)','m(z)','2(110)']}
3279    elif '-42m' in SytSym:
3280        ops = {'(x)':['-4(x)','2(y)','m(yz)'],'(y)':['-4(y)','2(z)','m(xz)'],'(z)':['-4(z)','2(x)','m(110)']}
3281    elif '-4' in SytSym:
3282        ops = {'(x)':['-4(x)',],'(y)':['-4(y)',],'(z)':['-4(z)',],}
3283    elif '4' in SytSym:
3284        ops = {'(x)':['4(x)',],'(y)':['4(y)',],'(z)':['4(z)',],}
3285
3286    elif '222' in SytSym:
3287        ops = {'':['2(x)','2(y)','2(z)'],
3288                   '(x)':['2(y)','2(z)','2(x)'],'(y)':['2(x)','2(z)','2(y)'],'(z)':['2(x)','2(y)','2(z)'],
3289                   '(100)':['2(z)','2(100)','2(120)',],'(010)':['2(z)','2(010)','2(210)',],
3290                   '(110)':['2(z)','2(110)','2(+-0)',],}
3291    elif 'mm2' in SytSym:
3292        ops = {'(x)':['m(y)','m(z)','2(x)'],'(y)':['m(x)','m(z)','2(y)'],'(z)':['m(x)','m(y)','2(z)'],
3293               '(xy)':['m(+-0)','m(z)','2(110)'],'(yz)':['m(0+-)','m(xz)','2(yz)'],     #not 2(xy)!
3294               '(xz)':['m(+0-)','m(y)','2(xz)'],'(z100)':['m(100)','m(120)','2(z)'],
3295               '(z010)':['m(010)','m(210)','2(z)'],'(z110)':['m(110)','m(+-0)','2(z)'],
3296               '(+-0)':[ 'm(110)','m(z)','2(+-0)'],'(d100)':['m(yz)','m(0+-)','2(xz)'],
3297               '(d010)':['m(xz)','m(+0-)','2(y)'],'(d001)':['m(110)','m(+-0)','2(z)'],
3298               '(210)':['m(z)','m(010)','2(210)'],'(120)':['m(z)','m(100)','2(120)'],
3299               '(100)':['m(z)','m(120)','2(100)',],'(010)':['m(z)','m(210)','2(010)',],
3300               '(110)':['m(z)','m(+-0)','2(110)',],}
3301    elif 'mmm' in SytSym:
3302        ops = {'':['m(x)','m(y)','m(z)'],
3303                   '(100)':['m(z)','m(100)','m(120)',],'(010)':['m(z)','m(010)','m(210)',],
3304                   '(110)':['m(z)','m(110)','m(+-0)',],
3305                   '(x)':['m(x)','m(y)','m(z)'],'(y)':['m(x)','m(y)','m(z)'],'(z)':['m(x)','m(y)','m(z)'],}
3306       
3307    elif '32' in SytSym:
3308        ops = {'(120)':['3','2(120)',],'(100)':['3','2(100)']}
3309    elif '23' in SytSym:
3310        ops = {'':['2(x)','3(111)']}
3311    elif 'm3' in SytSym:
3312        ops = {'(100)':['(+-0)',],'(+--)':[],'(-+-)':[],'(--+)':[]}
3313    elif '3m' in SytSym:
3314        ops = {'(111)':['3(111)','m(+-0)',],'(+--)':['3(+--)','m(0+-)',],
3315               '(-+-)':['3(-+-)','m(+0-)',],'(--+)':['3(--+)','m(+-0)',],
3316               '(100)':['3','m(100)'],'(120)':['3','m(210)',]}
3317   
3318    if SytSym.split('(')[0] in ['6/m','6mm','-6m2','622','-6','-3','-3m',]:     #not simple cases
3319        MagSytSym = SytSym
3320        if "-1'" in dupDir:
3321            if '-6' in SytSym:
3322                MagSytSym = MagSytSym.replace('-6',"-6'")
3323            elif '-3m' in SytSym:
3324                MagSytSym = MagSytSym.replace('-3m',"-3'm'")
3325            elif '-3' in SytSym:
3326                MagSytSym = MagSytSym.replace('-3',"-3'")
3327        elif '-6m2' in SytSym:
3328            if "m'(110)" in dupDir:
3329                MagSytSym = "-6m'2'("+SytSym.split('(')[1]
3330        elif '6/m' in SytSym:
3331            if "m'(z)" in dupDir:
3332                MagSytSym = "6'/m'"
3333        elif '6mm' in SytSym:
3334            if "m'(110)" in dupDir:
3335                MagSytSym = "6'm'm"
3336        return MagSytSym
3337    try:
3338        axis = '('+SytSym.split('(')[1]
3339    except IndexError:
3340        axis = ''
3341    MagSytSym = ''
3342    for m in ops[axis]:
3343        if m in dupDir:
3344            MagSytSym += m.split('(')[0]
3345        else:
3346            MagSytSym += m.split('(')[0]+"'"
3347        if '2/m' in SytSym and '2' in m:
3348            MagSytSym += '/'
3349        if '-3/m' in SytSym:
3350            MagSytSym = '-'+MagSytSym
3351       
3352    MagSytSym += axis
3353# some exceptions & special rules         
3354    if MagSytSym == "4'/m'm'm'": MagSytSym = "4/m'm'm'"
3355    return MagSytSym
3356   
3357#    if len(GenSym) == 3:
3358#        if SGSpin[1] < 0:
3359#            if 'mm2' in SytSym:
3360#                MagSytSym = "m'm'2"+'('+SplitSytSym[1]
3361#            else:   #bad rule for I41/a
3362#                MagSytSym = SplitSytSym[0]+"'"
3363#                if len(SplitSytSym) > 1:
3364#                    MagSytSym += '('+SplitSytSym[1]
3365#        else:
3366#            MagSytSym = SytSym
3367#        if len(SplitSytSym) >1:
3368#            if "-4'"+'('+SplitSytSym[1] in dupDir:
3369#                MagSytSym = MagSytSym.replace('-4',"-4'")
3370#            if "-6'"+'('+SplitSytSym[1] in dupDir:
3371#                MagSytSym = MagSytSym.replace('-6',"-6'")
3372#        return MagSytSym
3373#           
3374    return SytSym
3375
3376def UpdateSytSym(Phase):
3377    ''' Update site symmetry/site multiplicity after space group/BNS lattice change
3378    '''
3379    generalData = Phase['General']
3380    SGData = generalData['SGData']
3381    Atoms = Phase['Atoms']
3382    cx,ct,cs,cia = generalData['AtomPtrs']
3383    for atom in Atoms:
3384        XYZ = atom[cx:cx+3]
3385        sytsym,Mult = SytSym(XYZ,SGData)[:2]
3386        sytSym,Mul,Nop,dupDir = SytSym(XYZ,SGData)
3387        atom[cs] = sytsym
3388        if generalData['Type'] == 'magnetic':
3389            magSytSym = MagSytSym(sytSym,dupDir,SGData)
3390            atom[cs] = magSytSym
3391        atom[cs+1] = Mult
3392    return
3393   
3394def ElemPosition(SGData):
3395    ''' Under development.
3396    Object here is to return a list of symmetry element types and locations suitable
3397    for say drawing them.
3398    So far I have the element type... getting all possible locations without lookup may be impossible!
3399    '''
3400    Inv = SGData['SGInv']
3401    eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
3402    # get operators & expand if centrosymmetric
3403    Ops = SGData['SGOps']
3404    opM = np.array([op[0].T for op in Ops])
3405    opT = np.array([op[1] for op in Ops])
3406    if Inv:
3407        opM = np.concatenate((opM,-opM))
3408        opT = np.concatenate((opT,-opT))
3409    opMT = list(zip(opM,opT))
3410    for M,T in opMT[1:]:        #skip I
3411        Dt = int(nl.det(M))
3412        Tr = int(np.trace(M))
3413        Dt = -(Dt-1)//2
3414        Es = eleSym[Tr][Dt]
3415        if Dt:              #rotation-inversion
3416            I = np.eye(3)
3417            if Tr == 1:     #mirrors/glides
3418                if np.any(T):       #glide
3419                    M2 = np.inner(M,M)
3420                    MT = np.inner(M,T)+T
3421                    print ('glide',Es,MT)
3422                    print (M2)
3423                else:               #mirror
3424                    print ('mirror',Es,T)
3425                    print (I-M)
3426                X = [-1,-1,-1]
3427            elif Tr == -3:  # pure inversion
3428                X = np.inner(nl.inv(I-M),T)
3429                print ('inversion',Es,X)
3430            else:           #other rotation-inversion
3431                M2 = np.inner(M,M)
3432                MT = np.inner(M,T)+T
3433                print ('rot-inv',Es,MT)
3434                print (M2)
3435                X = [-1,-1,-1]
3436        else:               #rotations
3437            print ('rotation',Es)
3438            X = [-1,-1,-1]
3439        #SymElements.append([Es,X])
3440       
3441    return #SymElements
3442   
3443def ApplyStringOps(A,SGData,X,Uij=[]):
3444    'Needs a doc string'
3445    SGOps = SGData['SGOps']
3446    SGCen = SGData['SGCen']
3447    Ax = A.split('+')
3448    Ax[0] = int(Ax[0])
3449    iC = 1
3450    if Ax[0] < 0:
3451        iC = -1
3452    Ax[0] = abs(Ax[0])
3453    nA = Ax[0]%100-1
3454    cA = Ax[0]//100
3455    Cen = SGCen[cA]
3456    M,T = SGOps[nA]
3457    if len(Ax)>1:
3458        cellA = Ax[1].split(',')
3459        cellA = np.array([int(a) for a in cellA])
3460    else:
3461        cellA = np.zeros(3)
3462    newX = Cen+iC*(np.inner(M,X).T+T)+cellA
3463    if len(Uij):
3464        U = Uij2U(Uij)
3465        U = np.inner(M,np.inner(U,M).T)
3466        newUij = U2Uij(U)
3467        return [newX,newUij]
3468    else:
3469        return newX
3470       
3471def ApplyStringOpsMom(A,SGData,Mom):
3472    'Needs a doc string'
3473    SGOps = SGData['SGOps']
3474    Ax = A.split('+')
3475    Ax[0] = int(Ax[0])
3476    iAx = abs(Ax[0])
3477    nA = iAx%100-1
3478    nC = len(SGOps)*(iAx//100)
3479    NA = nA
3480    if Ax[0] < 0:
3481        NA += len(SGOps)
3482    M,T = SGOps[nA]
3483    if SGData['SGGray']:        #no nonzero moments for gray groups!
3484        newMom = [0.,0.,0.]
3485    else:
3486        newMom = np.inner(Mom,M).T*nl.det(M)*SGData['SpnFlp'][NA+nC]
3487#        print(len(SGOps),Ax[0],iAx,nC,nA,NA,SGData['SpnFlp'][NA],Mom,newMom)
3488    return newMom
3489       
3490def StringOpsProd(A,B,SGData):
3491    """
3492    Find A*B where A & B are in strings '-' + '100*c+n' + '+ijk'
3493    where '-' indicates inversion, c(>0) is the cell centering operator,
3494    n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3495    Should return resultant string - C. SGData - dictionary using entries:
3496
3497       *  'SGCen': cell centering vectors [0,0,0] at least
3498       *  'SGOps': symmetry operations as [M,T] so that M*x+T = x'
3499
3500    """
3501    SGOps = SGData['SGOps']
3502    SGCen = SGData['SGCen']
3503    #1st split out the cell translation part & work on the operator parts
3504    Ax = A.split('+'); Bx = B.split('+')
3505    Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
3506    iC = 0
3507    if Ax[0]*Bx[0] < 0:
3508        iC = 1
3509    Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
3510    nA = Ax[0]%100-1;  nB = Bx[0]%100-1
3511    cA = Ax[0]//100;  cB = Bx[0]//100
3512    Cen = (SGCen[cA]+SGCen[cB])%1.0
3513    cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
3514    Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
3515    Mc = np.inner(Ma,Mb.T)
3516#    print Ma,Mb,Mc
3517    Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
3518#    print Ta,Tb,Tc
3519#    print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
3520    nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
3521    #now the cell translation part
3522    if len(Ax)>1:
3523        cellA = Ax[1].split(',')
3524        cellA = [int(a) for a in cellA]
3525    else:
3526        cellA = [0,0,0]
3527    if len(Bx)>1:
3528        cellB = Bx[1].split(',')
3529        cellB = [int(b) for b in cellB]
3530    else:
3531        cellB = [0,0,0]
3532    cellC = np.add(cellA,cellB)
3533    C = str(((nC+1)+(100*cC))*(1-2*iC))+'+'+ \
3534        str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
3535    return C
3536           
3537def U2Uij(U):
3538    #returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
3539    return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
3540   
3541def Uij2U(Uij):
3542    #returns the thermal motion tensor U from Uij as numpy array
3543    return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
3544
3545def StandardizeSpcName(spcgroup):
3546    '''Accept a spacegroup name where spaces may have not been used
3547    in the names according to the GSAS convention (spaces between symmetry
3548    for each axis) and return the space group name as used in GSAS
3549    '''
3550    rspc = spcgroup.replace(' ','').upper()
3551    # deal with rhombohedral and hexagonal setting designations
3552    rhomb = ''
3553    if rspc[-1:] == 'R':
3554        rspc = rspc[:-1]
3555        rhomb = ' R'
3556    gray = ''
3557    if "1'" in rspc:
3558        gray = " 1'"
3559        rspc = rspc.replace("1'",'')
3560    rspc = rspc.replace("'",'')
3561    if rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
3562        rspc = rspc[:-1]
3563    if rspc[1:3] in ['M3','N3','A3','D3']:      #fix cubic old style
3564        rspc.replace('3','-3')
3565    bns = -1
3566    try:
3567        bns = rspc.index('_')
3568        rspc = rspc.replace(rspc[bns:bns+2],'')
3569    except ValueError:
3570        pass
3571    # look for a match in the spacegroup lists
3572    for i in spglist.values():
3573        for spc in i:
3574            if rspc == spc.replace(' ','').upper():
3575                return spc+gray+rhomb
3576    # how about the post-2002 orthorhombic names?
3577    if rspc in sgequiv_2002_orthorhombic:
3578        return sgequiv_2002_orthorhombic[rspc]+gray
3579    else:
3580    # not found
3581        return ''
3582   
3583def SpaceGroupNumber(spcgroup):
3584    SGNo = -1
3585    SpcGp = StandardizeSpcName(spcgroup)
3586    if not SpcGp:
3587        return SGNo
3588    try:
3589        SGNo = spgbyNum.index(SpcGp)
3590    except ValueError:
3591        pass
3592    return SGNo
3593
3594
3595spgbyNum = []
3596'''Space groups indexed by number'''
3597spgbyNum = [None,
3598        'P 1','P -1',                                                   #1-2
3599        'P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
3600        'C 2/m','P 2/c','P 21/c','C 2/c',                               #3-15
3601        'P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21',
3602        'C 2 2 21','C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
3603        'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
3604        'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
3605        'C m m 2','C m c 21','C c c 2',
3606        'A m m 2','A b m 2','A m a 2','A b a 2',
3607        'F m m 2','F d d 2','I m m 2','I b a 2','I m a 2',
3608        'P m m m','P n n n','P c c m','P b a n',
3609        'P m m a','P n n a','P m n a','P c c a','P b a m','P c c n',
3610        'P b c m','P n n m','P m m n','P b c n','P b c a','P n m a',
3611        'C m c m','C m c a','C m m m','C c c m','C m m a','C c c a',
3612        'F m m m', 'F d d d',
3613        'I m m m','I b a m','I b c a','I m m a',                        #16-74
3614        'P 4','P 41','P 42','P 43',
3615        'I 4','I 41',
3616        'P -4','I -4','P 4/m','P 42/m','P 4/n','P 42/n',
3617        'I 4/m','I 41/a',
3618        'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3619        'P 42 21 2','P 43 2 2','P 43 21 2',
3620        'I 4 2 2','I 41 2 2',
3621        'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c','P 4 n c',
3622        'P 42 m c','P 42 b c',
3623        'I 4 m m','I 4 c m','I 41 m d','I 41 c d',
3624        'P -4 2 m','P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2',
3625        'P -4 c 2','P -4 b 2','P -4 n 2',
3626        'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d',
3627        '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',
3628        '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',
3629        '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',
3630        'P 42/n c m',
3631        'I 4/m m m','I 4/m c m','I 41/a m d','I 41/a c d',
3632        'P 3','P 31','P 32','R 3','P -3','R -3',
3633        'P 3 1 2','P 3 2 1','P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1',
3634        'R 3 2',
3635        'P 3 m 1','P 3 1 m','P 3 c 1','P 3 1 c',
3636        'R 3 m','R 3 c',
3637        'P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1',
3638        'R -3 m','R -3 c',                                               #75-167
3639        'P 6','P 61',
3640        'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3641        'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3642        'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3643        '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
3644        'P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m 3','P n 3',
3645        'F m -3','F d -3','I m -3',
3646        'P a -3','I a -3','P 4 3 2','P 42 3 2','F 4 3 2','F 41 3 2',
3647        'I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2','P -4 3 m',
3648        'F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
3649        'P m -3 m','P n -3 n','P m -3 n','P n -3 m',
3650        'F m -3 m','F m -3 c','F d -3 m','F d -3 c',