source: trunk/GSASIIspc.py @ 3589

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

fixes to Bilbao mag symm stuff - now gets nonstd. sp.grp. correct with right site symm rules

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