source: trunk/GSASIIspc.py @ 3380

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

Add several text blocks to tutorial listing & switch format to display them.
fix bug involving new restraint layout & an old gpx file.
add GetLittleGrpOps? routine to G2spc

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