source: trunk/GSASIIspc.py @ 3396

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

set refList=False for a couple of calls to GetHistogramPhaseData? where reflection list shouldn't be rebuilt
fix OnDrawDist? to report Sym Ops correctly
Move atom labels over 4 spaces to be clear of atom ball

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 153.8 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-22 20:51:10 +0000 (Tue, 22 May 2018) $
12# $Author: vondreele $
13# $Revision: 3396 $
14# $URL: trunk/GSASIIspc.py $
15# $Id: GSASIIspc.py 3396 2018-05-22 20:51:10Z 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: 3396 $")
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 SGProd(OpA,OpB):
1730    '''
1731    Form space group operator product. OpA & OpB are [M,V] pairs;
1732        both must be of same dimension (3 or 4). Returns [M,V] pair
1733    '''
1734    A,U = OpA
1735    B,V = OpB
1736    M = np.inner(B,A.T)
1737    W = np.inner(B,U)+V
1738    return M,W
1739       
1740def GetLittleGrpOps(SGData,vec):
1741    ''' Find rotation part of operators that leave vec unchanged
1742   
1743    :param SGData: space group data structure as defined in SpcGroup above.
1744    :param vec: a numpy array of fractional vector coordinates
1745    :returns: Little - list of operators [M,T] that form the little gropu
1746    '''
1747    Little = []
1748    Ops = SGData['SGOps'][:]
1749    if SGData['SGInv']:
1750        Ops += [[-M,-T] for [M,T] in Ops]
1751    for [M,T] in Ops:
1752        tvec = np.inner(M,vec)%1.
1753        if np.allclose(tvec,vec%1.):
1754            Little.append([M,T])
1755    return Little
1756       
1757def MoveToUnitCell(xyz):
1758    '''
1759    Translates a set of coordinates so that all values are >=0 and < 1
1760
1761    :param xyz: a list or numpy array of fractional coordinates
1762    :returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
1763    '''
1764    XYZ = (np.array(xyz)+10.)%1.
1765    cell = np.asarray(np.rint(XYZ-xyz),dtype=np.int32)
1766    return XYZ,cell
1767       
1768def Opposite(XYZ,toler=0.0002):
1769    '''
1770    Gives opposite corner, edge or face of unit cell for position within tolerance.
1771        Result may be just outside the cell within tolerance
1772
1773    :param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
1774    :param toler: unit cell fraction tolerance making opposite
1775    :returns:
1776        XYZ: dict of opposite positions; key=unit cell & always contains XYZ
1777    '''
1778    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]]
1779    TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
1780    perm = TB*perm3
1781    cperm = ['%d,%d,%d'%(i,j,k) for i,j,k in perm]
1782    D = dict(zip(cperm,perm))
1783    new = {}
1784    for key in D:
1785        new[key] = np.array(D[key])+np.array(XYZ)
1786    return new
1787       
1788def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
1789    '''
1790    Generates the equivalent positions for a specified coordinate and space group
1791
1792    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
1793    :param SGData: from :func:`SpcGroup`
1794    :param All: True return all equivalent positions including duplicates;
1795      False return only unique positions
1796    :param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
1797    :param Move: True move generated atom positions to be inside cell
1798      False do not move atoms       
1799    :return: [[XYZEquiv],Idup,[UijEquiv],spnflp]
1800
1801      *  [XYZEquiv] is list of equivalent positions (XYZ is first entry)
1802      *  Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
1803      * is centering operator number (1-4) and - is for inversion
1804        Cell = unit cell translations needed to put new positions inside cell
1805        [UijEquiv] - equivalent Uij; absent if no Uij given
1806      * +1/-1 for spin inversion of operator - empty if not magnetic
1807       
1808    '''
1809    XYZEquiv = []
1810    UijEquiv = []
1811    Idup = []
1812    Cell = []
1813    inv = int(SGData['SGInv']+1)
1814    icen = SGData['SGCen']
1815    if SGData['SGFixed']:
1816        inv = 1
1817    SpnFlp = SGData.get('SpnFlp',[])
1818    spnflp = []
1819    X = np.array(XYZ)
1820    mj = 0
1821    if Move:
1822        X = MoveToUnitCell(X)[0]
1823    for ic,cen in enumerate(icen):
1824        C = np.array(cen)
1825        for invers in range(inv):
1826            for io,[M,T] in enumerate(SGData['SGOps']):
1827                idup = ((io+1)+100*ic)*(1-2*invers)
1828                XT = np.inner(M,X)+T
1829                if len(Uij):
1830                    U = Uij2U(Uij)
1831                    U = np.inner(M,np.inner(U,M).T)
1832                    newUij = U2Uij(U)
1833                if invers:
1834                    XT = -XT
1835                XT += C
1836                cell = np.zeros(3,dtype=np.int32)
1837                cellj = np.zeros(3,dtype=np.int32)
1838                if Move:
1839                    newX,cellj = MoveToUnitCell(XT)
1840                else:
1841                    newX = XT
1842                cell += cellj
1843                if All:
1844                    if np.allclose(newX,X,atol=0.0002):
1845                        idup = False
1846                else:
1847                    if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
1848                        idup = False
1849                if All or idup:
1850                    XYZEquiv.append(newX)
1851                    Idup.append(idup)
1852                    Cell.append(cell)
1853                    if len(Uij):
1854                        UijEquiv.append(newUij)
1855                    if len(SpnFlp):
1856                        spnflp.append(SpnFlp[mj])
1857                    else:
1858                        spnflp.append(1)
1859                mj += 1
1860    if len(Uij):
1861        return zip(XYZEquiv,UijEquiv,Idup,Cell,spnflp)
1862    else:
1863        return zip(XYZEquiv,Idup,Cell,spnflp)
1864       
1865def GenHKL(HKL,SGData):
1866    ''' Generates all equivlent reflections including Friedel pairs
1867    :param HKL:  [h,k,l] must be integral values
1868    :param SGData: space group data obtained from SpcGroup
1869    :returns: array Uniq: equivalent reflections
1870    '''
1871   
1872    Ops = SGData['SGOps']
1873    OpM = np.array([op[0] for op in Ops])
1874    Uniq = np.inner(OpM,HKL)
1875    Uniq = list(Uniq)+list(-1*Uniq)
1876    return np.array(Uniq)
1877
1878def GenHKLf(HKL,SGData):
1879    '''
1880    Uses old GSAS Fortran routine genhkl.for
1881
1882    :param HKL:  [h,k,l] must be integral values for genhkl.for to work
1883    :param SGData: space group data obtained from SpcGroup
1884    :returns: iabsnt,mulp,Uniq,phi
1885
1886     *   iabsnt = True if reflection is forbidden by symmetry
1887     *   mulp = reflection multiplicity including Friedel pairs
1888     *   Uniq = numpy array of equivalent hkl in descending order of h,k,l
1889     *   phi = phase offset for each equivalent h,k,l
1890
1891    '''
1892    hklf = list(HKL)+[0,]       #could be numpy array!
1893    Ops = SGData['SGOps']
1894    OpM = np.array([op[0] for op in Ops],order='F')
1895    OpT = np.array([op[1] for op in Ops])
1896    Cen = np.array([cen for cen in SGData['SGCen']],order='F')
1897   
1898    import pyspg
1899    Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
1900    h,k,l,f = Uniq
1901    Uniq=np.array(list(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq])))
1902    phi = f[:Nuniq]
1903    return iabsnt,mulp,Uniq,phi
1904
1905#def MagHKLchk(HKL,SGData):
1906#    SpnFlp = SGData['SpnFlp']
1907#    print(HKL)
1908#    Uniq = GenHKL(HKL,SGData)
1909   
1910def checkSSLaue(HKL,SGData,SSGData):
1911    #Laue check here - Toss HKL if outside unique Laue part
1912    h,k,l,m = HKL
1913    if SGData['SGLaue'] == '2/m':
1914        if SGData['SGUniq'] == 'a':
1915            if 'a' in SSGData['modSymb'] and h == 0 and m < 0:
1916                return False
1917            elif 'b' in SSGData['modSymb'] and k == 0 and l ==0 and m < 0:
1918                return False
1919            else:
1920                return True
1921        elif SGData['SGUniq'] == 'b':
1922            if 'b' in SSGData['modSymb'] and k == 0 and m < 0:
1923                return False
1924            elif 'a' in SSGData['modSymb'] and h == 0 and l ==0 and m < 0:
1925                return False
1926            else:
1927                return True
1928        elif SGData['SGUniq'] == 'c':
1929            if 'g' in SSGData['modSymb'] and l == 0 and m < 0:
1930                return False
1931            elif 'a' in SSGData['modSymb'] and h == 0 and k ==0 and m < 0:
1932                return False
1933            else:
1934                return True
1935    elif SGData['SGLaue'] == 'mmm':
1936        if 'a' in SSGData['modSymb']:
1937            if h == 0 and m < 0:
1938                return False
1939            else:
1940                return True
1941        elif 'b' in SSGData['modSymb']:
1942            if k == 0 and m < 0:
1943                return False
1944            else:
1945                return True
1946        elif 'g' in SSGData['modSymb']:
1947            if l == 0 and m < 0:
1948                return False
1949            else:
1950                return True
1951    else:   #tetragonal, trigonal, hexagonal (& triclinic?)
1952        if l == 0 and m < 0:
1953            return False
1954        else:
1955            return True
1956       
1957   
1958def checkSSextc(HKL,SSGData):
1959    Ops = SSGData['SSGOps']
1960    OpM = np.array([op[0] for op in Ops])
1961    OpT = np.array([op[1] for op in Ops])
1962    HKLS = np.array([HKL,-HKL])     #Freidel's Law
1963    DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,4))
1964    PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
1965    for dhkl,phkl in zip(DHKL,PHKL)[1:]:    #skip identity
1966        if dhkl.any():
1967            continue
1968        else:
1969            if phkl%1.:
1970                return False
1971    return True
1972   
1973################################################################################
1974#### Site symmetry tables
1975################################################################################
1976     
1977OprName = {
1978    '-6643':       ['-1',1],'6479' :    ['2(z)',2],'-6479':     ['m(z)',3],
1979    '6481' :     ['m(y)',4],'-6481':    ['2(y)',5],'6641' :     ['m(x)',6],
1980    '-6641':     ['2(x)',7],'6591' :  ['m(+-0)',8],'-6591':   ['2(+-0)',9],
1981    '6531' :  ['m(110)',10],'-6531': ['2(110)',11],'6537' :    ['4(z)',12],
1982    '-6537':   ['-4(z)',13],'975'  : ['3(111)',14],'6456' :       ['3',15],
1983    '-489' :  ['3(+--)',16],'483'  : ['3(-+-)',17],'-969' :  ['3(--+)',18],
1984    '819'  :  ['m(+0-)',19],'-819' : ['2(+0-)',20],'2431' :  ['m(0+-)',21],
1985    '-2431':  ['2(0+-)',22],'-657' :  ['m(xz)',23],'657'  :   ['2(xz)',24],
1986    '1943' :   ['-4(x)',25],'-1943':   ['4(x)',26],'-2429':   ['m(yz)',27],
1987    '2429' :   ['2(yz)',28],'639'  :  ['-4(y)',29],'-639' :    ['4(y)',30],
1988    '-6484':   ['2(010)',4],'6484' :  ['m(010)',5],'-6668':   ['2(100)',6],
1989    '6668' :   ['m(100)',7],'-6454': ['2(120)',18],'6454' :  ['m(120)',19],
1990    '-6638':  ['2(210)',20],'6638' : ['m(210)',21],   #search in SytSym ends at m(210)
1991    '2223' : ['3(+++)2',39],
1992    '6538' :   ['6(z)1',40],'-2169':['3(--+)2',41],'2151' : ['3(+--)2',42],
1993    '2205' :['-3(-+-)2',43],'-2205':[' (-+-)2',44],'489'  :['-3(+--)1',45],
1994    '801'  :   ['4(y)1',46],'1945' :  ['4(x)3',47],'-6585': ['-4(z)3 ',48],
1995    '6585' :   ['4(z)3',49],'6584' :  ['3(z)2',50],'6666' :  ['6(z)5 ',51],
1996    '6643' :       ['1',52],'-801' : ['-4(y)1',53],'-1945': ['-4(x)3 ',54],
1997    '-6666':  ['-6(z)5',55],'-6538': ['-6(z)1',56],'-2223':['-3(+++)2',57],
1998    '-975' :['-3(+++)1',58],'-6456': ['-3(z)1',59],'-483' :['-3(-+-)1',60],
1999    '969'  :['-3(--+)1',61],'-6584': ['-3(z)2',62],'2169' :['-3(--+)2',63],
2000    '-2151':['-3(+--)2',64],   }                               
2001
2002KNsym = {
2003    '0'         :'    1   ','1'         :'   -1   ','64'        :'    2(x)','32'        :'    m(x)',
2004    '97'        :'  2/m(x)','16'        :'    2(y)','8'         :'    m(y)','25'        :'  2/m(y)',
2005    '2'         :'    2(z)','4'         :'    m(z)','7'         :'  2/m(z)','134217728' :'   2(yz)',
2006    '67108864'  :'   m(yz)','201326593' :' 2/m(yz)','2097152'   :'  2(0+-)','1048576'   :'  m(0+-)',
2007    '3145729'   :'2/m(0+-)','8388608'   :'   2(xz)','4194304'   :'   m(xz)','12582913'  :' 2/m(xz)',
2008    '524288'    :'  2(+0-)','262144'    :'  m(+0-)','796433'    :'2/m(+0-)','1024'      :'   2(xy)',
2009    '512'       :'   m(xy)','1537'      :' 2/m(xy)','256'       :'  2(+-0)','128'       :'  m(+-0)',
2010    '385'       :'2/m(+-0)','76'        :'  mm2(x)','52'        :'  mm2(y)','42'        :'  mm2(z)',
2011    '135266336' :' mm2(yz)','69206048'  :'mm2(0+-)','8650760'   :' mm2(xz)','4718600'   :'mm2(+0-)',
2012    '1156'      :' mm2(xy)','772'       :'mm2(+-0)','82'        :'  222   ','136314944' :'  222(x)',
2013    '8912912'   :'  222(y)','1282'      :'  222(z)','127'       :'  mmm   ','204472417' :'  mmm(x)',
2014    '13369369'  :'  mmm(y)','1927'      :'  mmm(z)','33554496'  :'  4(x)','16777280'  :' -4(x)',
2015    '50331745'  :'4/m(x)'  ,'169869394' :'422(x)','84934738'  :'-42m(x)','101711948' :'4mm(x)',
2016    '254804095' :'4/mmm(x)','536870928 ':'  4(y)','268435472' :' -4(y)','805306393' :'4/m(y)',
2017    '545783890' :'422(y)','272891986' :'-42m(y)','541327412' :'4mm(y)','818675839' :'4/mmm(y)',
2018    '2050'      :'  4(z)','4098'      :' -4(z)','6151'      :'4/m(z)','3410'      :'422(z)',
2019    '4818'      :'-42m(z)','2730'      :'4mm(z)','8191'      :'4/mmm(z)','8192'      :'  3(111)',
2020    '8193'      :' -3(111)','2629888'   :' 32(111)','1319040'   :' 3m(111)','3940737'   :'-3m(111)',
2021    '32768'     :'  3(+--)','32769'     :' -3(+--)','10519552'  :' 32(+--)','5276160'   :' 3m(+--)',
2022    '15762945'  :'-3m(+--)','65536'     :'  3(-+-)','65537'     :' -3(-+-)','134808576' :' 32(-+-)',
2023    '67437056'  :' 3m(-+-)','202180097' :'-3m(-+-)','131072'    :'  3(--+)','131073'    :' -3(--+)',
2024    '142737664' :' 32(--+)','71434368'  :' 3m(--+)','214040961' :'-3m(--+)','237650'    :'   23   ',
2025    '237695'    :'   m3   ','715894098' :'   432  ','358068946' :'  -43m  ','1073725439':'   m3m  ',
2026    '68157504'  :' mm2(d100)','4456464'   :' mm2(d010)','642'       :' mm2(d001)','153092172' :'-4m2(x)',
2027    '277348404' :'-4m2(y)','5418'      :'-4m2(z)','1075726335':'  6/mmm ','1074414420':'-6m2(100)',
2028    '1075070124':'-6m2(120)','1075069650':'   6mm  ','1074414890':'   622  ','1073758215':'   6/m  ',
2029    '1073758212':'   -6   ','1073758210':'    6   ','1073759865':'-3m(100)','1075724673':'-3m(120)',
2030    '1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
2031    '1073758209':'   -3   ','1073758208':'    3   ','1074135143':'mmm(100)','1075314719':'mmm(010)',
2032    '1073743751':'mmm(110)','1074004034':' mm2(z100)','1074790418':' mm2(z010)','1073742466':' mm2(z110)',
2033    '1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
2034    '1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
2035    '1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
2036    '1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
2037    '1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
2038    '1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
2039    '1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
2040    '1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':'   -1   ',
2041    }
2042
2043NXUPQsym = {
2044    '1'        :(28,29,28,28),'-1'       :( 1,29,28, 0),'2(x)'     :(12,18,12,25),'m(x)'     :(25,18,12,25),
2045    '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),
2046    '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),
2047    '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),
2048    '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),
2049    '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),
2050    '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),
2051    '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),
2052    'mm2(yz)'  :(10,13, 0,-1),'mm2(0+-)' :(11,13, 0,-1),'mm2(xz)'  :( 8,12, 0,-1),'mm2(+0-)' :( 9,12, 0,-1),
2053    'mm2(xy)'  :( 6,11, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222'      :( 1,10, 0,-1),'222(x)'   :( 1,13, 0,-1),
2054    '222(y)'   :( 1,12, 0,-1),'222(z)'   :( 1,11, 0,-1),'mmm'      :( 1,10, 0,-1),'mmm(x)'   :( 1,13, 0,-1),
2055    'mmm(y)'   :( 1,12, 0,-1),'mmm(z)'   :( 1,11, 0,-1),'4(x)'     :(12, 4,12, 0),'-4(x)'    :( 1, 4,12, 0),
2056    '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),
2057    '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),
2058    '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,),
2059    '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),
2060    '-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),
2061    '-3(111)'  :( 1, 5, 2, 0),'32(111)'  :( 1, 5, 0, 2),'3m(111)'  :( 2, 5, 0, 2),'-3m(111)' :( 1, 5, 0,-1),
2062    '3(+--)'   :( 5, 8, 5, 0),'-3(+--)'  :( 1, 8, 5, 0),'32(+--)'  :( 1, 8, 0, 5),'3m(+--)'  :( 5, 8, 0, 5),
2063    '-3m(+--)' :( 1, 8, 0,-1),'3(-+-)'   :( 4, 7, 4, 0),'-3(-+-)'  :( 1, 7, 4, 0),'32(-+-)'  :( 1, 7, 0, 4),
2064    '3m(-+-)'  :( 4, 7, 0, 4),'-3m(-+-)' :( 1, 7, 0,-1),'3(--+)'   :( 3, 6, 3, 0),'-3(--+)'  :( 1, 6, 3, 0),
2065    '32(--+)'  :( 1, 6, 0, 3),'3m(--+)'  :( 3, 6, 0, 3),'-3m(--+)' :( 1, 6, 0,-1),'23'       :( 1, 1, 0, 0),
2066    'm3'       :( 1, 1, 0, 0),'432'      :( 1, 1, 0, 0),'-43m'     :( 1, 1, 0, 0),'m3m'      :( 1, 1, 0, 0),
2067    'mm2(d100)':(12,13, 0,-1),'mm2(d010)':(13,12, 0,-1),'mm2(d001)':(14,11, 0,-1),'-4m2(x)'  :( 1, 4, 0,-1),
2068    '-4m2(y)'  :( 1, 3, 0,-1),'-4m2(z)'  :( 1, 2, 0,-1),'6/mmm'    :( 1, 9, 0,-1),'-6m2(100)':( 1, 9, 0,-1),
2069    '-6m2(120)':( 1, 9, 0,-1),'6mm'      :(14, 9, 0,-1),'622'      :( 1, 9, 0,-1),'6/m'      :( 1, 9,14,-1),
2070    '-6'       :( 1, 9,14, 0),'6'        :(14, 9,14, 0),'-3m(100)' :( 1, 9, 0,-1),'-3m(120)' :( 1, 9, 0,-1),
2071    '3m(100)'  :(14, 9, 0,14),'3m(120)'  :(14, 9, 0,14),'32(100)'  :( 1, 9, 0,14),'32(120)'  :( 1, 9, 0,14),
2072    '-3'       :( 1, 9,14, 0),'3'        :(14, 9,14, 0),'mmm(100)' :( 1,14, 0,-1),'mmm(010)' :( 1,15, 0,-1),
2073    'mmm(110)' :( 1,11, 0,-1),'mm2(z100)':(14,14, 0,-1),'mm2(z010)':(14,15, 0,-1),'mm2(z110)':(14,11, 0,-1),
2074    'mm2(100)' :(12,14, 0,-1),'mm2(010)' :(13,15, 0,-1),'mm2(110)' :( 6,11, 0,-1),'mm2(120)' :(15,14, 0,-1),
2075    'mm2(210)' :(16,15, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222(100)' :( 1,14, 0,-1),'222(010)' :( 1,15, 0,-1),
2076    '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),
2077    '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),
2078    'm(001)'   :(23,16,14,23),'m(100)'   :(26,25,12,26),'m(010)'   :(27,28,13,27),'m(110)'   :(18,19, 6,18),
2079    'm(120)'   :(24,27,15,24),'m(210)'   :(25,26,16,25),'m(+-0)'   :(17,20, 7,17),'2(001)'   :(14,16,14,23),
2080    '2(100)'   :(12,25,12,26),'2(010)'   :(13,28,13,27),'2(110)'   :( 6,19, 6,18),'2(120)'   :(15,27,15,24),
2081    '2(210)'   :(16,26,16,25),'2(+-0)'   :( 7,20, 7,17),'-1'       :( 1,29,28, 0)
2082    }
2083       
2084CSxinel = [[],      # 0th empty - indices are Fortran style
2085    [[0,0,0],[ 0.0, 0.0, 0.0]],      #1  0  0  0
2086    [[1,1,1],[ 1.0, 1.0, 1.0]],      #2  X  X  X
2087    [[1,1,1],[ 1.0, 1.0,-1.0]],      #3  X  X -X
2088    [[1,1,1],[ 1.0,-1.0, 1.0]],      #4  X -X  X
2089    [[1,1,1],[ 1.0,-1.0,-1.0]],      #5 -X  X  X
2090    [[1,1,0],[ 1.0, 1.0, 0.0]],      #6  X  X  0
2091    [[1,1,0],[ 1.0,-1.0, 0.0]],      #7  X -X  0
2092    [[1,0,1],[ 1.0, 0.0, 1.0]],      #8  X  0  X
2093    [[1,0,1],[ 1.0, 0.0,-1.0]],      #9  X  0 -X
2094    [[0,1,1],[ 0.0, 1.0, 1.0]],      #10  0  Y  Y
2095    [[0,1,1],[ 0.0, 1.0,-1.0]],      #11 0  Y -Y
2096    [[1,0,0],[ 1.0, 0.0, 0.0]],      #12  X  0  0
2097    [[0,1,0],[ 0.0, 1.0, 0.0]],      #13  0  Y  0
2098    [[0,0,1],[ 0.0, 0.0, 1.0]],      #14  0  0  Z
2099    [[1,1,0],[ 1.0, 2.0, 0.0]],      #15  X 2X  0
2100    [[1,1,0],[ 2.0, 1.0, 0.0]],      #16 2X  X  0
2101    [[1,1,2],[ 1.0, 1.0, 1.0]],      #17  X  X  Z
2102    [[1,1,2],[ 1.0,-1.0, 1.0]],      #18  X -X  Z
2103    [[1,2,1],[ 1.0, 1.0, 1.0]],      #19  X  Y  X
2104    [[1,2,1],[ 1.0, 1.0,-1.0]],      #20  X  Y -X
2105    [[1,2,2],[ 1.0, 1.0, 1.0]],      #21  X  Y  Y
2106    [[1,2,2],[ 1.0, 1.0,-1.0]],      #22  X  Y -Y
2107    [[1,2,0],[ 1.0, 1.0, 0.0]],      #23  X  Y  0
2108    [[1,0,2],[ 1.0, 0.0, 1.0]],      #24  X  0  Z
2109    [[0,1,2],[ 0.0, 1.0, 1.0]],      #25  0  Y  Z
2110    [[1,1,2],[ 1.0, 2.0, 1.0]],      #26  X 2X  Z
2111    [[1,1,2],[ 2.0, 1.0, 1.0]],      #27 2X  X  Z
2112    [[1,2,3],[ 1.0, 1.0, 1.0]],      #28  X  Y  Z
2113    ]
2114
2115CSuinel = [[],      # 0th empty - indices are Fortran style
2116    [[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
2117    [[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
2118    [[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
2119    [[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
2120    [[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
2121    [[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
2122    [[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
2123    [[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
2124    [[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
2125    [[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
2126    [[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
2127    [[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
2128    [[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
2129    [[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
2130    [[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
2131    [[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
2132    [[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
2133    [[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
2134    [[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
2135    [[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
2136    [[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
2137    [[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
2138    [[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
2139    [[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
2140    [[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
2141    [[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
2142    [[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
2143    [[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
2144    [[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
2145    ]
2146   
2147################################################################################
2148#### Site symmetry routines
2149################################################################################
2150   
2151def GetOprPtrName(key):
2152    'Needs a doc string'
2153    oprName = OprName[key][0]
2154    return oprName.replace('(','').replace(')','')
2155
2156def GetOprPtrNumber(key):
2157    'Needs a doc string'
2158    return OprName[key][1]
2159
2160def GetOprName(key):
2161    'Needs a doc string'
2162    return OprName[key][0]
2163
2164def GetKNsym(key):
2165    'Needs a doc string'
2166    try:
2167        return KNsym[key].strip()
2168    except KeyError:
2169        return 'sp'
2170
2171def GetNXUPQsym(siteSym):
2172    '''       
2173    The codes XUPQ are for lookup of symmetry constraints for position(X), thermal parm(U) & magnetic moments (P & Q)
2174    '''
2175    return NXUPQsym[siteSym]
2176
2177def GetCSxinel(siteSym): 
2178    "returns Xyz terms, multipliers, GUI flags"
2179    indx = GetNXUPQsym(siteSym.strip())
2180    return CSxinel[indx[0]]
2181   
2182def GetCSuinel(siteSym):
2183    "returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
2184    indx = GetNXUPQsym(siteSym.strip())
2185    return CSuinel[indx[1]]
2186   
2187def GetCSpqinel(SpnFlp,dupDir): 
2188    "returns Mxyz terms, multipliers, GUI flags"
2189    CSI = [[1,2,3],[1.0,1.0,1.0]]
2190    for sopr in dupDir:
2191#        print (sopr,dupDir[sopr])
2192        opr = sopr.replace("'",'')
2193        indx = GetNXUPQsym(opr)
2194        if SpnFlp[dupDir[sopr]] > 0:
2195            csi = CSxinel[indx[2]]  #P
2196        else:
2197            csi = CSxinel[indx[3]]  #Q
2198#        print(opr,indx,csi,CSI)
2199        if not len(csi):
2200            return [[0,0,0],[0.,0.,0.]]
2201        for kcs in [0,1,2]:
2202            if csi[0][kcs] == 0 and CSI[0][kcs] != 0:
2203                jcs = CSI[0][kcs]
2204                for ics in [0,1,2]:
2205                    if CSI[0][ics] == jcs:
2206                        CSI[0][ics] = 0
2207                        CSI[1][ics] = 0.
2208                    elif CSI[0][ics] > jcs:
2209                        CSI[0][ics] = CSI[0][ics]-1
2210            elif (CSI[0][kcs] == csi[0][kcs]) and (CSI[1][kcs] != csi[1][kcs]):
2211                CSI[1][kcs] = csi[1][kcs]
2212            elif CSI[0][kcs] >= csi[0][kcs]:
2213                CSI[0][kcs] = min(CSI[0][kcs],csi[0][kcs])
2214                if CSI[1][kcs] != csi[1][kcs]:
2215                    if CSI[1][kcs] == 1.:
2216                        CSI[1][kcs] = csi[1][kcs]
2217#        print(CSI)
2218    return CSI
2219   
2220def getTauT(tau,sop,ssop,XYZ,wave=np.zeros(3)):
2221    phase = np.sum(XYZ*wave)
2222    ssopinv = nl.inv(ssop[0])
2223    mst = ssopinv[3][:3]
2224    epsinv = ssopinv[3][3]
2225    sdet = nl.det(sop[0])
2226    ssdet = nl.det(ssop[0])
2227    dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
2228    dT = 1.0
2229    if np.any(dtau%.5):
2230        sumdtau = np.sum(dtau%.5)
2231        dT = 0.
2232        if np.abs(sumdtau-.5) > 1.e-4:
2233            dT = np.tan(np.pi*sumdtau)
2234    tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]+phase)
2235    return sdet,ssdet,dtau,dT,tauT
2236   
2237def OpsfromStringOps(A,SGData,SSGData):
2238    SGOps = SGData['SGOps']
2239    SSGOps = SSGData['SSGOps']
2240    Ax = A.split('+')
2241    Ax[0] = int(Ax[0])
2242    iC = 1
2243    if Ax[0] < 0:
2244        iC = -1
2245    Ax[0] = abs(Ax[0])
2246    nA = Ax[0]%100-1
2247    nC = Ax[0]//100
2248    unit = [0,0,0]
2249    if len(Ax) > 1:
2250        unit = eval('['+Ax[1]+']')
2251    return SGOps[nA],SSGOps[nA],iC,SGData['SGCen'][nC],unit
2252   
2253def GetSSfxuinel(waveType,Stype,nH,XYZ,SGData,SSGData,debug=False):
2254   
2255    def orderParms(CSI):
2256        parms = [0,]
2257        for csi in CSI:
2258            for i in [0,1,2]:
2259                if csi[i] not in parms:
2260                    parms.append(csi[i])
2261        for csi in CSI:
2262            for i in [0,1,2]:
2263                csi[i] = parms.index(csi[i])
2264        return CSI
2265       
2266    def fracCrenel(tau,Toff,Twid):
2267        Tau = (tau-Toff[:,np.newaxis])%1.
2268        A = np.where(Tau<Twid[:,np.newaxis],1.,0.)
2269        return A
2270       
2271    def fracFourier(tau,nH,fsin,fcos):
2272        SA = np.sin(2.*nH*np.pi*tau)
2273        CB = np.cos(2.*nH*np.pi*tau)
2274        A = SA[np.newaxis,np.newaxis,:]*fsin[:,:,np.newaxis]
2275        B = CB[np.newaxis,np.newaxis,:]*fcos[:,:,np.newaxis]
2276        return A+B
2277       
2278    def posFourier(tau,nH,psin,pcos):
2279        SA = np.sin(2*nH*np.pi*tau)
2280        CB = np.cos(2*nH*np.pi*tau)
2281        A = SA[np.newaxis,np.newaxis,:]*psin[:,:,np.newaxis]
2282        B = CB[np.newaxis,np.newaxis,:]*pcos[:,:,np.newaxis]
2283        return A+B   
2284
2285    def posSawtooth(tau,Toff,slopes):
2286        Tau = (tau-Toff)%1.
2287        A = slopes[:,np.newaxis]*Tau
2288        return A
2289   
2290    def posZigZag(tau,Tmm,XYZmax):
2291        DT = Tmm[1]-Tmm[0]
2292        slopeUp = 2.*XYZmax/DT
2293        slopeDn = 2.*XYZmax/(1.-DT)
2294        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])
2295        return A
2296
2297    def posBlock(tau,Tmm,XYZmax):
2298        A = np.array([np.where(Tmm[0] < t <= Tmm[1],XYZmax,-XYZmax) for t in tau])
2299        return A
2300       
2301    def DoFrac():
2302        dF = None
2303        dFTP = None
2304        if siteSym == '1':
2305            CSI = [[1,0],[2,0]],2*[[1.,0.],]
2306        elif siteSym == '-1':
2307            CSI = [[1,0],[0,0]],2*[[1.,0.],]
2308        else:
2309            delt2 = np.eye(2)*0.001
2310            FSC = np.ones(2,dtype='i')
2311            CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
2312            if 'Crenel' in waveType:
2313                dF = np.zeros_like(tau)
2314            else:
2315                dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2316            dFT = np.zeros_like(dF)
2317            dFTP = []
2318            for i in SdIndx:
2319                sop = Sop[i]
2320                ssop = SSop[i]           
2321                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2322                fsc = np.ones(2,dtype='i')
2323                if 'Crenel' in waveType:
2324                    dFT = np.zeros_like(tau)
2325                    fsc = [1,1]
2326                else:   #Fourier
2327                    dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
2328                    dFT = nl.det(sop[0])*dFT
2329                    dFT = dFT[:,np.argsort(tauT)]
2330                    dFT[0] *= ssdet
2331                    dFT[1] *= sdet
2332                    dFTP.append(dFT)
2333               
2334                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2335                        fsc = [1,1]
2336                        if dT:
2337                            CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
2338                        else:
2339                            CSI = [[[1,0],[0,0]],[[1.,0.],[0.,0.]]]
2340                        FSC = np.zeros(2,dtype='i')
2341                        return CSI,dF,dFTP
2342                    else:
2343                        for i in range(2):
2344                            if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
2345                                fsc[i] = 1
2346                            else:
2347                                fsc[i] = 0
2348                        FSC &= fsc
2349                        if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc)
2350            n = -1
2351            for i,F in enumerate(FSC):
2352                if F:
2353                    n += 1
2354                    CSI[0][i] = n+1
2355                    CSI[1][i] = 1.0
2356           
2357        return CSI,dF,dFTP
2358       
2359    def DoXYZ():
2360        dX,dXTP = None,None
2361        if siteSym == '1':
2362            CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2363        elif siteSym == '-1':
2364            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.],]
2365        else:
2366            delt4 = np.ones(4)*0.001
2367            delt5 = np.ones(5)*0.001
2368            delt6 = np.eye(6)*0.001
2369            if 'Fourier' in waveType:
2370                dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2371                  #3x6x12 modulated position array (X,Spos,tau)& force positive
2372                CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2373            elif waveType == 'Sawtooth':
2374                dX = posSawtooth(tau,delt4[0],delt4[1:])
2375                CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0]]),
2376                    np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2377            elif waveType in ['ZigZag','Block']:
2378                if waveType == 'ZigZag':
2379                    dX = posZigZag(tau,delt5[:2],delt5[2:])
2380                else:
2381                    dX = posBlock(tau,delt5[:2],delt5[2:])
2382                CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2383                    np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2384            XSC = np.ones(6,dtype='i')
2385            dXTP = []
2386            for i in SdIndx:
2387                sop = Sop[i]
2388                ssop = SSop[i]
2389                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2390                xsc = np.ones(6,dtype='i')
2391                if 'Fourier' in waveType:
2392                    dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2393                elif waveType == 'Sawtooth':
2394                    dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2395                elif waveType == 'ZigZag':
2396                    dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2397                elif waveType == 'Block':
2398                    dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2399                dXT = np.inner(sop[0],dXT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2400                dXT = np.swapaxes(dXT,1,2)      # back to array(3x6x49)
2401                dXT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2402                dXTP.append(dXT)
2403                if waveType == 'Fourier':
2404                    for i in range(3):
2405                        if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2406                            xsc[i] = 0
2407                        if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2408                            xsc[i+3] = 0
2409                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2410                        xsc[3:6] = 0
2411                        CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2412                            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2413                        if dT:
2414                            if '(x)' in siteSym:
2415                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2416                                if 'm' in siteSym and len(SdIndx) == 1:
2417                                    CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2418                            elif '(y)' in siteSym:
2419                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2420                                if 'm' in siteSym and len(SdIndx) == 1:
2421                                    CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2422                            elif '(z)' in siteSym:
2423                                CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2424                                if 'm' in siteSym and len(SdIndx) == 1:
2425                                    CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2426                        else:
2427                            CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2428                    if '4/mmm' in laue:
2429                        if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2430                            if '(xy)' in siteSym:
2431                                CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2432                                if dT:
2433                                    CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2434                                else:
2435                                    CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2436                        if '(xy)' in siteSym or '(+-0)' in siteSym:
2437                            mul = 1
2438                            if '(+-0)' in siteSym:
2439                                mul = -1
2440                            if np.allclose(dX[0,0,:],dXT[1,0,:]):
2441                                CSI[0][3:5] = [[11,0,0],[11,0,0]]
2442                                CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2443                                xsc[3:5] = 0
2444                            if np.allclose(dX[0,3,:],dXT[0,4,:]):
2445                                CSI[0][:2] = [[12,0,0],[12,0,0]]
2446                                CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2447                                xsc[:2] = 0
2448                XSC &= xsc
2449                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2450            if waveType == 'Fourier':
2451                n = -1
2452                if debug: print (XSC)
2453                for i,X in enumerate(XSC):
2454                    if X:
2455                        n += 1
2456                        CSI[0][i][0] = n+1
2457                        CSI[1][i][0] = 1.0
2458           
2459        return list(CSI),dX,dXTP
2460       
2461    def DoUij():
2462        dU,dUTP = None,None
2463        if siteSym == '1':
2464            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2465                [7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2466        elif siteSym == '-1':
2467            CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],   \
2468                6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2469        else:
2470            tau = np.linspace(0,1,49,True)
2471            delt12 = np.eye(12)*0.0001
2472            dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
2473            CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2474            USC = np.ones(12,dtype='i')
2475            dUTP = []
2476            dtau = 0.
2477            for i in SdIndx:
2478                sop = Sop[i]
2479                ssop = SSop[i]
2480                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2481                usc = np.ones(12,dtype='i')
2482                dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
2483                dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3)    #convert dUT to 12x49x3x3
2484                dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2485                dUT = np.array(U2Uij(dUijT))    #convert to 6x12x49
2486                dUT = dUT[:,:,np.argsort(tauT)]
2487                dUT[:,:6,:] *=(ssdet*sdet)
2488                dUTP.append(dUT)
2489                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2490                    if dT:
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                        [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2495                    else:
2496                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2497                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2498                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2499                        [0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2500                    if 'mm2(x)' in siteSym and dT:
2501                        CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2502                        USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2503                    elif '(xy)' in siteSym and dT:
2504                        CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2505                            [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2506                        CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2507                        USC = [1,1,1,1,1,1,1,1,1,1,1,1]                             
2508                    elif '(x)' in siteSym and dT:
2509                        CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2510                    elif '(y)' in siteSym and dT:
2511                        CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2512                    elif '(z)' in siteSym and dT:
2513                        CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2514                    for i in range(6):
2515                        if not USC[i]:
2516                            CSI[0][i] = [0,0,0]
2517                            CSI[1][i] = [0.,0.,0.]
2518                            CSI[0][i+6] = [0,0,0]
2519                            CSI[1][i+6] = [0.,0.,0.]
2520                else:                       
2521                    for i in range(6):
2522                        if not np.allclose(dU[i,i,:],dUT[i,i,:]):  #sin part
2523                            usc[i] = 0
2524                        if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]):   #cos part
2525                            usc[i+6] = 0
2526                    if np.any(dUT[1,0,:]):
2527                        if '4/m' in siteSym:
2528                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2529                            if ssop[1][3]:
2530                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2531                                usc[9] = 1
2532                            else:
2533                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2534                                usc[9] = 0
2535                        elif '4' in siteSym:
2536                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2537                            CSI[0][:2] = [[11,0,0],[11,0,0]]
2538                            if ssop[1][3]:
2539                                CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2540                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2541                                usc[2] = 0
2542                                usc[8] = 0
2543                                usc[3] = 1
2544                                usc[9] = 1
2545                            else:
2546                                CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2547                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2548                                usc[2] = 1
2549                                usc[8] = 1
2550                                usc[3] = 0               
2551                                usc[9] = 0
2552                        elif 'xy' in siteSym or '+-0' in siteSym:
2553                            if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2554                                CSI[0][4:6] = [[12,0,0],[12,0,0]]
2555                                CSI[0][6:8] = [[11,0,0],[11,0,0]]
2556                                CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2557                                CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2558                                usc[4:6] = 0
2559                                usc[6:8] = 0
2560                           
2561                    if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2562                USC &= usc
2563            if debug: print (USC)
2564            if not np.any(dtau%.5):
2565                n = -1
2566                for i,U in enumerate(USC):
2567                    if U:
2568                        n += 1
2569                        CSI[0][i][0] = n+1
2570                        CSI[1][i][0] = 1.0
2571   
2572        return list(CSI),dU,dUTP
2573   
2574    def DoMag():
2575        CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2576        if siteSym == '1':
2577            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2578        elif siteSym in ['-1','mmm','2/m(x)','2/m(y)','2/m(z)','4/mmm001']:
2579            CSI = 3*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3*[[0.,0.,0.],]+3*[[1.,0.,0.],]
2580        else:
2581            delt6 = np.eye(6)*0.001
2582            dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,np.newaxis,np.newaxis]
2583              #3x6x12 modulated moment array (M,Spos,tau)& force positive
2584            CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2585            MSC = np.ones(6,dtype='i')
2586            dMTP = []
2587            for i in SdIndx:
2588                sop = Sop[i]
2589                ssop = SSop[i]
2590                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2591                msc = np.ones(6,dtype='i')
2592                dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2593                dMT = np.inner(sop[0],dMT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2594                dMT = np.swapaxes(dMT,1,2)      # back to array(3x6x49)
2595                dMT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2596                dMTP.append(dMT)
2597                for i in range(3):
2598                    if not np.allclose(dM[i,i,:],dMT[i,i,:]):
2599                        msc[i] = 0
2600                    if not np.allclose(dM[i,i+3,:],dMT[i,i+3,:]):
2601                        msc[i+3] = 0
2602                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2603                    msc[3:6] = 0
2604                    CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2605                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2606                    if dT:
2607                        if '(x)' in siteSym:
2608                            CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2609                            if 'm' in siteSym and len(SdIndx) == 1:
2610                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2611                        elif '(y)' in siteSym:
2612                            CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2613                            if 'm' in siteSym and len(SdIndx) == 1:
2614                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2615                        elif '(z)' in siteSym:
2616                            CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2617                            if 'm' in siteSym and len(SdIndx) == 1:
2618                                CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2619                    else:
2620                        CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2621                if '4/mmm' in laue:
2622                    if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2623                        if '(xy)' in siteSym:
2624                            CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2625                            if dT:
2626                                CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2627                            else:
2628                                CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2629                    if '(xy)' in siteSym or '(+-0)' in siteSym:
2630                        mul = 1
2631                        if '(+-0)' in siteSym:
2632                            mul = -1
2633                        if np.allclose(dM[0,0,:],dMT[1,0,:]):
2634                            CSI[0][3:5] = [[11,0,0],[11,0,0]]
2635                            CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2636                            msc[3:5] = 0
2637                        if np.allclose(dM[0,3,:],dMT[0,4,:]):
2638                            CSI[0][:2] = [[12,0,0],[12,0,0]]
2639                            CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2640                            msc[:2] = 0
2641                MSC &= msc
2642                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2643            n = -1
2644            if debug: print (MSC)
2645            for i,M in enumerate(MSC):
2646                if M:
2647                    n += 1
2648                    CSI[0][i][0] = n+1
2649                    CSI[1][i][0] = 1.0
2650
2651        return CSI,None,None
2652       
2653    if debug: print ('super space group: '+SSGData['SSpGrp'])
2654    xyz = np.array(XYZ)%1.
2655    SGOps = copy.deepcopy(SGData['SGOps'])
2656    laue = SGData['SGLaue']
2657    siteSym = SytSym(XYZ,SGData)[0].strip()
2658    if debug: print ('siteSym: '+siteSym)
2659    SSGOps = copy.deepcopy(SSGData['SSGOps'])
2660    #expand ops to include inversions if any
2661    if SGData['SGInv']:
2662        for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
2663            SGOps.append([-op[0],-op[1]%1.])
2664            SSGOps.append([-sop[0],-sop[1]%1.])
2665    #build set of sym ops around special position       
2666    SSop = []
2667    Sop = []
2668    Sdtau = []
2669    for iop,Op in enumerate(SGOps):         
2670        nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
2671        if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
2672            SSop.append(SSGOps[iop])
2673            Sop.append(SGOps[iop])
2674            ssopinv = nl.inv(SSGOps[iop][0])
2675            mst = ssopinv[3][:3]
2676            epsinv = ssopinv[3][3]
2677            Sdtau.append(np.sum(mst*(XYZ-SGOps[iop][1])-epsinv*SSGOps[iop][1][3]))
2678    SdIndx = np.argsort(np.array(Sdtau))     # just to do in sensible order
2679    if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
2680    #setup displacement arrays
2681    tau = np.linspace(-1,1,49,True)
2682    #make modulation arrays - one parameter at a time
2683    if Stype == 'Sfrac':
2684        CSI,dF,dFTP = DoFrac()
2685    elif Stype == 'Spos':
2686        CSI,dF,dFTP = DoXYZ()
2687        CSI[0] = orderParms(CSI[0])
2688    elif Stype == 'Sadp':
2689        CSI,dF,dFTP = DoUij()
2690        CSI[0] = orderParms(CSI[0]) 
2691    elif Stype == 'Smag':
2692        CSI,dF,dFTP = DoMag()
2693       
2694    if debug:
2695        return CSI,dF,dFTP
2696    else:
2697        return CSI
2698   
2699def MustrainNames(SGData):
2700    'Needs a doc string'
2701    laue = SGData['SGLaue']
2702    uniq = SGData['SGUniq']
2703    if laue in ['m3','m3m']:
2704        return ['S400','S220']
2705    elif laue in ['6/m','6/mmm','3m1']:
2706        return ['S400','S004','S202']
2707    elif laue in ['31m','3']:
2708        return ['S400','S004','S202','S211']
2709    elif laue in ['3R','3mR']:
2710        return ['S400','S220','S310','S211']
2711    elif laue in ['4/m','4/mmm']:
2712        return ['S400','S004','S220','S022']
2713    elif laue in ['mmm']:
2714        return ['S400','S040','S004','S220','S202','S022']
2715    elif laue in ['2/m']:
2716        SHKL = ['S400','S040','S004','S220','S202','S022']
2717        if uniq == 'a':
2718            SHKL += ['S013','S031','S211']
2719        elif uniq == 'b':
2720            SHKL += ['S301','S103','S121']
2721        elif uniq == 'c':
2722            SHKL += ['S130','S310','S112']
2723        return SHKL
2724    else:
2725        SHKL = ['S400','S040','S004','S220','S202','S022']
2726        SHKL += ['S310','S103','S031','S130','S301','S013']
2727        SHKL += ['S211','S121','S112']
2728        return SHKL
2729       
2730def HStrainVals(HSvals,SGData):
2731    laue = SGData['SGLaue']
2732    uniq = SGData['SGUniq']
2733    DIJ = np.zeros(6)
2734    if laue in ['m3','m3m']:
2735        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
2736    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2737        DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
2738    elif laue in ['3R','3mR']:
2739        DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
2740    elif laue in ['4/m','4/mmm']:
2741        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
2742    elif laue in ['mmm']:
2743        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2744    elif laue in ['2/m']:
2745        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2746        if uniq == 'a':
2747            DIJ[5] = HSvals[3]
2748        elif uniq == 'b':
2749            DIJ[4] = HSvals[3]
2750        elif uniq == 'c':
2751            DIJ[3] = HSvals[3]
2752    else:
2753        DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
2754    return DIJ
2755
2756def HStrainNames(SGData):
2757    'Needs a doc string'
2758    laue = SGData['SGLaue']
2759    uniq = SGData['SGUniq']
2760    if laue in ['m3','m3m']:
2761        return ['D11','eA']         #add cubic strain term
2762    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2763        return ['D11','D33']
2764    elif laue in ['3R','3mR']:
2765        return ['D11','D12']
2766    elif laue in ['4/m','4/mmm']:
2767        return ['D11','D33']
2768    elif laue in ['mmm']:
2769        return ['D11','D22','D33']
2770    elif laue in ['2/m']:
2771        Dij = ['D11','D22','D33']
2772        if uniq == 'a':
2773            Dij += ['D23']
2774        elif uniq == 'b':
2775            Dij += ['D13']
2776        elif uniq == 'c':
2777            Dij += ['D12']
2778        return Dij
2779    else:
2780        Dij = ['D11','D22','D33','D12','D13','D23']
2781        return Dij
2782   
2783def MustrainCoeff(HKL,SGData):
2784    'Needs a doc string'
2785    #NB: order of terms is the same as returned by MustrainNames
2786    laue = SGData['SGLaue']
2787    uniq = SGData['SGUniq']
2788    h,k,l = HKL
2789    Strm = []
2790    if laue in ['m3','m3m']:
2791        Strm.append(h**4+k**4+l**4)
2792        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
2793    elif laue in ['6/m','6/mmm','3m1']:
2794        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
2795        Strm.append(l**4)
2796        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
2797    elif laue in ['31m','3']:
2798        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
2799        Strm.append(l**4)
2800        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
2801        Strm.append(4.0*h*k*l*(h+k))
2802    elif laue in ['3R','3mR']:
2803        Strm.append(h**4+k**4+l**4)
2804        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
2805        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
2806        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
2807    elif laue in ['4/m','4/mmm']:
2808        Strm.append(h**4+k**4)
2809        Strm.append(l**4)
2810        Strm.append(3.0*(h*k)**2)
2811        Strm.append(3.0*((h*l)**2+(k*l)**2))
2812    elif laue in ['mmm']:
2813        Strm.append(h**4)
2814        Strm.append(k**4)
2815        Strm.append(l**4)
2816        Strm.append(3.0*(h*k)**2)
2817        Strm.append(3.0*(h*l)**2)
2818        Strm.append(3.0*(k*l)**2)
2819    elif laue in ['2/m']:
2820        Strm.append(h**4)
2821        Strm.append(k**4)
2822        Strm.append(l**4)
2823        Strm.append(3.0*(h*k)**2)
2824        Strm.append(3.0*(h*l)**2)
2825        Strm.append(3.0*(k*l)**2)
2826        if uniq == 'a':
2827            Strm.append(2.0*k*l**3)
2828            Strm.append(2.0*l*k**3)
2829            Strm.append(4.0*k*l*h**2)
2830        elif uniq == 'b':
2831            Strm.append(2.0*l*h**3)
2832            Strm.append(2.0*h*l**3)
2833            Strm.append(4.0*h*l*k**2)
2834        elif uniq == 'c':
2835            Strm.append(2.0*h*k**3)
2836            Strm.append(2.0*k*h**3)
2837            Strm.append(4.0*h*k*l**2)
2838    else:
2839        Strm.append(h**4)
2840        Strm.append(k**4)
2841        Strm.append(l**4)
2842        Strm.append(3.0*(h*k)**2)
2843        Strm.append(3.0*(h*l)**2)
2844        Strm.append(3.0*(k*l)**2)
2845        Strm.append(2.0*k*h**3)
2846        Strm.append(2.0*h*l**3)
2847        Strm.append(2.0*l*k**3)
2848        Strm.append(2.0*h*k**3)
2849        Strm.append(2.0*l*h**3)
2850        Strm.append(2.0*k*l**3)
2851        Strm.append(4.0*k*l*h**2)
2852        Strm.append(4.0*h*l*k**2)
2853        Strm.append(4.0*k*h*l**2)
2854    return Strm
2855
2856def MuShklMean(SGData,Amat,Shkl):
2857   
2858    def genMustrain(xyz,Shkl):
2859        uvw = np.inner(Amat.T,xyz)
2860        Strm = np.array(MustrainCoeff(uvw,SGData))
2861        Sum = np.sum(np.multiply(Shkl,Strm))
2862        Sum = np.where(Sum > 0.01,Sum,0.01)
2863        Sum = np.sqrt(Sum)
2864        return Sum*xyz
2865       
2866    PHI = np.linspace(0.,360.,30,True)
2867    PSI = np.linspace(0.,180.,30,True)
2868    X = np.outer(npcosd(PHI),npsind(PSI))
2869    Y = np.outer(npsind(PHI),npsind(PSI))
2870    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
2871    XYZ = np.dstack((X,Y,Z))
2872    XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
2873    return np.sqrt(np.sum(XYZ**2)/900.)
2874   
2875def Muiso2Shkl(muiso,SGData,cell):
2876    "this is to convert isotropic mustrain to generalized Shkls"
2877    import GSASIIlattice as G2lat
2878    A = G2lat.cell2AB(cell)[0]
2879   
2880    def minMus(Shkl,muiso,H,SGData,A):
2881        U = np.inner(A.T,H)
2882        S = np.array(MustrainCoeff(U,SGData))
2883        Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,np.newaxis]),axis=0))
2884        rad = np.sqrt(np.sum((Sum[:,np.newaxis]*H)**2,axis=1))
2885        return (muiso-rad)**2
2886       
2887    laue = SGData['SGLaue']
2888    PHI = np.linspace(0.,360.,60,True)
2889    PSI = np.linspace(0.,180.,60,True)
2890    X = np.outer(npsind(PHI),npsind(PSI))
2891    Y = np.outer(npcosd(PHI),npsind(PSI))
2892    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
2893    HKL = np.dstack((X,Y,Z))
2894    if laue in ['m3','m3m']:
2895        S0 = [1000.,1000.]
2896    elif laue in ['6/m','6/mmm','3m1']:
2897        S0 = [1000.,1000.,1000.]
2898    elif laue in ['31m','3']:
2899        S0 = [1000.,1000.,1000.,1000.]
2900    elif laue in ['3R','3mR']:
2901        S0 = [1000.,1000.,1000.,1000.]
2902    elif laue in ['4/m','4/mmm']:
2903        S0 = [1000.,1000.,1000.,1000.]
2904    elif laue in ['mmm']:
2905        S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
2906    elif laue in ['2/m']:
2907        S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
2908    else:
2909        S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000., 
2910            1000.,1000.,0.,0.,0.]
2911    S0 = np.array(S0)
2912    HKL = np.reshape(HKL,(-1,3))
2913    result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
2914    return result[0]
2915       
2916def PackRot(SGOps):
2917    IRT = []
2918    for ops in SGOps:
2919        M = ops[0]
2920        irt = 0
2921        for j in range(2,-1,-1):
2922            for k in range(2,-1,-1):
2923                irt *= 3
2924                irt += M[k][j]
2925        IRT.append(int(irt))
2926    return IRT
2927       
2928def SytSym(XYZ,SGData):
2929    '''
2930    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
2931
2932    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
2933    :param SGData: from SpcGroup
2934    :Returns: a four element tuple:
2935
2936     * The 1st element is a code for the site symmetry (see GetKNsym)
2937     * The 2nd element is the site multiplicity
2938     * Ndup number of overlapping operators
2939     * dupDir Dict - dictionary of overlapping operators
2940
2941    '''
2942    Mult = 1
2943    Isym = 0
2944    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
2945        Isym = 1073741824
2946    Jdup = 0
2947    Ndup = 0
2948    dupDir = {}
2949    inv = SGData['SGInv']+1
2950    icen = SGData['SGCen']
2951    if SGData['SGFixed']:       #already in list of operators
2952        inv = 1
2953    Xeqv = list(GenAtom(XYZ,SGData,True))
2954#    for xeqv in Xeqv:   print(xeqv)
2955    IRT = PackRot(SGData['SGOps'])
2956    L = -1
2957    for ic,cen in enumerate(icen):
2958        for invers in range(int(inv)):
2959            for io,ops in enumerate(SGData['SGOps']):
2960                irtx = (1-2*invers)*IRT[io]
2961                L += 1
2962                if not Xeqv[L][1]:
2963                    Ndup = io
2964                    Jdup += 1
2965                    jx = GetOprPtrNumber(str(irtx))   #[KN table no,op name,KNsym ptr]
2966                    if jx < 39:
2967                        px = GetOprName(str(irtx))
2968                        if Xeqv[L][-1] < 0:
2969                            if '(' in px:
2970                                px = px.split('(')
2971                                px[0] += "'"
2972                                px = '('.join(px)
2973                            else:   
2974                                px += "'"
2975                        dupDir[px] = L
2976                        Isym += 2**(jx-1)
2977    if Isym == 1073741824: Isym = 0
2978    Mult = len(SGData['SGOps'])*len(SGData['SGCen'])*inv//Jdup
2979         
2980    return GetKNsym(str(Isym)),Mult,Ndup,dupDir
2981   
2982def MagSytSym(SytSym,dupDir,SGData):
2983    '''
2984    site sym operations: 1,-1,2,3,-3,4,-4,6,-6,m need to be marked if spin inversion
2985    '''
2986    SGData['GenSym'],SGData['GenFlg'] = GetGenSym(SGData)[:2]
2987#    print('SGPtGrp',SGData['SGPtGrp'],'SytSym',SytSym,'MagSpGrp',SGData['MagSpGrp'])
2988#    print('dupDir',dupDir)
2989    SplitSytSym = SytSym.split('(')
2990    if SGData['SGGray']:
2991        return SytSym+"1'"
2992    if SytSym == '1':       #genersl position
2993        return SytSym
2994    if SplitSytSym[0] == SGData['SGPtGrp']:     #simple cases
2995        try:
2996            MagSytSym = SGData['MagSpGrp'].split()[1]
2997        except IndexError:
2998            MagSytSym = SGData['MagSpGrp'][1:].strip("1'")
2999        if len(SplitSytSym) > 1:
3000            MagSytSym += '('+SplitSytSym[1]
3001        return MagSytSym
3002    if len(dupDir) == 1:
3003        return dupDir.keys()[0]
3004   
3005   
3006    if '2/m' in SytSym:         #done I think; last 2wo might be not needed
3007        ops = {'(x)':['2(x)','m(x)'],'(y)':['2(y)','m(y)'],'(z)':['2(z)','m(z)'],
3008               '(100)':['2(100)','m(100)'],'(010)':['2(010)','m(010)'],'(001)':['2(001)','m(001)'],
3009               '(120)':['2(120)','m(120)'],'(210)':['2(210)','m(210)'],'(+-0)':['2(+-0)','m(+-0)'],
3010               '(110)':['2(110)','m(110)']}
3011   
3012    elif '4/mmm' in SytSym:
3013        ops = {'(x)':['4(x)','m(x)','m(y)','m(0+-)'],   #m(0+-) for cubic m3m?
3014               '(y)':['4(y)','m(y)','m(z)','m(+0-)'],   #m(+0-)
3015               '(z)':['4(z)','m(z)','m(x)','m(+-0)']}   #m(+-0)
3016    elif '4mm' in SytSym:
3017        ops = {'(x)':['4(x)','m(y)','m(yz)'],'(y)':['4(y)','m(z)','m(xz)'],'(z)':['4(z)','m(x)','m(110)']}
3018    elif '422' in SytSym:
3019        ops = {'(x)':['4(x)','2(y)','2(yz)'],'(y)':['4(y)','2(z)','2(xz)'],'(z)':['4(z)','2(x)','2(110)']}
3020    elif '-4m2' in SytSym:
3021        ops = {'(x)':['-4(x)','m(x)','2(yz)'],'(y)':['-4(y)','m(y)','2(xz)'],'(z)':['-4(z)','m(z)','2(110)']}
3022    elif '-42m' in SytSym:
3023        ops = {'(x)':['-4(x)','2(y)','m(yz)'],'(y)':['-4(y)','2(z)','m(xz)'],'(z)':['-4(z)','2(x)','m(110)']}
3024    elif '-4' in SytSym:
3025        ops = {'(x)':['-4(x)',],'(y)':['-4(y)',],'(z)':['-4(z)',],}
3026    elif '4' in SytSym:
3027        ops = {'(x)':['4(x)',],'(y)':['4(y)',],'(z)':['4(z)',],}
3028
3029    elif '222' in SytSym:
3030        ops = {'':['2(x)','2(y)','2(z)'],
3031                   '(x)':['2(y)','2(z)','2(x)'],'(y)':['2(x)','2(z)','2(y)'],'(z)':['2(x)','2(y)','2(z)'],
3032                   '(100)':['2(z)','2(100)','2(120)',],'(010)':['2(z)','2(010)','2(210)',],
3033                   '(110)':['2(z)','2(110)','2(+-0)',],}
3034    elif 'mm2' in SytSym:
3035        ops = {'(x)':['m(y)','m(z)','2(x)'],'(y)':['m(x)','m(z)','2(y)'],'(z)':['m(x)','m(y)','2(z)'],
3036               '(xy)':['m(+-0)','m(z)','2(110)'],'(yz)':['m(0+-)','m(xz)','2(yz)'],     #not 2(xy)!
3037               '(xz)':['m(+0-)','m(y)','2(xz)'],'(z100)':['m(100)','m(120)','2(z)'],
3038               '(z010)':['m(010)','m(210)','2(z)'],'(z110)':['m(110)','m(+-0)','2(z)'],
3039               '(+-0)':[ 'm(110)','m(z)','2(+-0)'],'(d100)':['m(yz)','m(0+-)','2(xz)'],
3040               '(d010)':['m(xz)','m(+0-)','2(y)'],'(d001)':['m(110)','m(+-0)','2(z)'],
3041               '(210)':['m(z)','m(010)','2(210)'],
3042               '(100)':['m(z)','m(120)','2(100)',],'(010)':['m(z)','m(210)','2(010)',],
3043               '(110)':['m(z)','m(+-0)','2(110)',],}
3044    elif 'mmm' in SytSym:
3045        ops = {'':['m(x)','m(y)','m(z)'],
3046                   '(100)':['m(z)','m(100)','m(120)',],'(010)':['m(z)','m(010)','m(210)',],
3047                   '(110)':['m(z)','m(110)','m(+-0)',],
3048                   '(x)':['m(x)','m(y)','m(z)'],'(y)':['m(x)','m(y)','m(z)'],'(z)':['m(x)','m(y)','m(z)'],}
3049       
3050    elif '32' in SytSym:
3051        ops = {'(120)':['3','2(120)',],'(100)':['3','2(100)']}
3052    elif '23' in SytSym:
3053        ops = {'':['2(x)','3(111)']}
3054    elif 'm3' in SytSym:
3055        ops = {'(100)':['(+-0)',],'(+--)':[],'(-+-)':[],'(--+)':[]}
3056    elif '3m' in SytSym:
3057        ops = {'(111)':['3(111)','m(+-0)',],'(+--)':['3(+--)','m(0+-)',],
3058               '(-+-)':['3(-+-)','m(+0-)',],'(--+)':['3(--+)','m(+-0)',],
3059               '(100)':['3','m(100)'],'(120)':['3','m(210)',]}
3060   
3061    if SytSym.split('(')[0] in ['6/m','6mm','-6m2','622','-6','-3','-3m',]:     #not simple cases
3062        MagSytSym = SytSym
3063        if "-1'" in dupDir:
3064            if '-6' in SytSym:
3065                MagSytSym = MagSytSym.replace('-6',"-6'")
3066            elif '-3m' in SytSym:
3067                MagSytSym = MagSytSym.replace('-3m',"-3'm'")
3068            elif '-3' in SytSym:
3069                MagSytSym = MagSytSym.replace('-3',"-3'")
3070        elif '-6m2' in SytSym:
3071            if "m'(110)" in dupDir:
3072                MagSytSym = "-6m'2'("+SytSym.split('(')[1]
3073        elif '6/m' in SytSym:
3074            if "m'(z)" in dupDir:
3075                MagSytSym = "6'/m'"
3076        elif '6mm' in SytSym:
3077            if "m'(110)" in dupDir:
3078                MagSytSym = "6'm'm"
3079        return MagSytSym
3080    try:
3081        axis = '('+SytSym.split('(')[1]
3082    except IndexError:
3083        axis = ''
3084    MagSytSym = ''
3085    for m in ops[axis]:
3086        if m in dupDir:
3087            MagSytSym += m.split('(')[0]
3088        else:
3089            MagSytSym += m.split('(')[0]+"'"
3090        if '2/m' in SytSym and '2' in m:
3091            MagSytSym += '/'
3092        if '-3/m' in SytSym:
3093            MagSytSym = '-'+MagSytSym
3094       
3095    MagSytSym += axis
3096# some exceptions & special rules         
3097    if MagSytSym == "4'/m'm'm'": MagSytSym = "4/m'm'm'"
3098    return MagSytSym
3099   
3100#    if len(GenSym) == 3:
3101#        if SGSpin[1] < 0:
3102#            if 'mm2' in SytSym:
3103#                MagSytSym = "m'm'2"+'('+SplitSytSym[1]
3104#            else:   #bad rule for I41/a
3105#                MagSytSym = SplitSytSym[0]+"'"
3106#                if len(SplitSytSym) > 1:
3107#                    MagSytSym += '('+SplitSytSym[1]
3108#        else:
3109#            MagSytSym = SytSym
3110#        if len(SplitSytSym) >1:
3111#            if "-4'"+'('+SplitSytSym[1] in dupDir:
3112#                MagSytSym = MagSytSym.replace('-4',"-4'")
3113#            if "-6'"+'('+SplitSytSym[1] in dupDir:
3114#                MagSytSym = MagSytSym.replace('-6',"-6'")
3115#        return MagSytSym
3116#           
3117    return SytSym
3118   
3119def ElemPosition(SGData):
3120    ''' Under development.
3121    Object here is to return a list of symmetry element types and locations suitable
3122    for say drawing them.
3123    So far I have the element type... getting all possible locations without lookup may be impossible!
3124    '''
3125    Inv = SGData['SGInv']
3126    eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
3127    # get operators & expand if centrosymmetric
3128    Ops = SGData['SGOps']
3129    opM = np.array([op[0].T for op in Ops])
3130    opT = np.array([op[1] for op in Ops])
3131    if Inv:
3132        opM = np.concatenate((opM,-opM))
3133        opT = np.concatenate((opT,-opT))
3134    opMT = list(zip(opM,opT))
3135    for M,T in opMT[1:]:        #skip I
3136        Dt = int(nl.det(M))
3137        Tr = int(np.trace(M))
3138        Dt = -(Dt-1)//2
3139        Es = eleSym[Tr][Dt]
3140        if Dt:              #rotation-inversion
3141            I = np.eye(3)
3142            if Tr == 1:     #mirrors/glides
3143                if np.any(T):       #glide
3144                    M2 = np.inner(M,M)
3145                    MT = np.inner(M,T)+T
3146                    print ('glide',Es,MT)
3147                    print (M2)
3148                else:               #mirror
3149                    print ('mirror',Es,T)
3150                    print (I-M)
3151                X = [-1,-1,-1]
3152            elif Tr == -3:  # pure inversion
3153                X = np.inner(nl.inv(I-M),T)
3154                print ('inversion',Es,X)
3155            else:           #other rotation-inversion
3156                M2 = np.inner(M,M)
3157                MT = np.inner(M,T)+T
3158                print ('rot-inv',Es,MT)
3159                print (M2)
3160                X = [-1,-1,-1]
3161        else:               #rotations
3162            print ('rotation',Es)
3163            X = [-1,-1,-1]
3164        #SymElements.append([Es,X])
3165       
3166    return #SymElements
3167   
3168def ApplyStringOps(A,SGData,X,Uij=[]):
3169    'Needs a doc string'
3170    SGOps = SGData['SGOps']
3171    SGCen = SGData['SGCen']
3172    Ax = A.split('+')
3173    Ax[0] = int(Ax[0])
3174    iC = 1
3175    if Ax[0] < 0:
3176        iC = -1
3177    Ax[0] = abs(Ax[0])
3178    nA = Ax[0]%100-1
3179    cA = Ax[0]//100
3180    Cen = SGCen[cA]
3181    M,T = SGOps[nA]
3182    if len(Ax)>1:
3183        cellA = Ax[1].split(',')
3184        cellA = np.array([int(a) for a in cellA])
3185    else:
3186        cellA = np.zeros(3)
3187    newX = Cen+iC*(np.inner(M,X).T+T)+cellA
3188    if len(Uij):
3189        U = Uij2U(Uij)
3190        U = np.inner(M,np.inner(U,M).T)
3191        newUij = U2Uij(U)
3192        return [newX,newUij]
3193    else:
3194        return newX
3195       
3196def ApplyStringOpsMom(A,SGData,Mom):
3197    'Needs a doc string'
3198    SGOps = SGData['SGOps']
3199    Ax = A.split('+')
3200    Ax[0] = int(Ax[0])
3201    iC = 1
3202    if Ax[0] < 0:
3203        iC = -1
3204    Ax[0] = abs(Ax[0])
3205    nA = Ax[0]%100-1
3206    M,T = SGOps[nA]
3207    if SGData['SGGray']:
3208        newMom = -(np.inner(Mom,M).T)*nl.det(M)
3209    else:
3210        newMom = -np.inner(Mom,M).T*SGData['MagMom'][nA-1]
3211    return newMom
3212       
3213def StringOpsProd(A,B,SGData):
3214    """
3215    Find A*B where A & B are in strings '-' + '100*c+n' + '+ijk'
3216    where '-' indicates inversion, c(>0) is the cell centering operator,
3217    n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3218    Should return resultant string - C. SGData - dictionary using entries:
3219
3220       *  'SGCen': cell centering vectors [0,0,0] at least
3221       *  'SGOps': symmetry operations as [M,T] so that M*x+T = x'
3222
3223    """
3224    SGOps = SGData['SGOps']
3225    SGCen = SGData['SGCen']
3226    #1st split out the cell translation part & work on the operator parts
3227    Ax = A.split('+'); Bx = B.split('+')
3228    Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
3229    iC = 0
3230    if Ax[0]*Bx[0] < 0:
3231        iC = 1
3232    Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
3233    nA = Ax[0]%100-1;  nB = Bx[0]%100-1
3234    cA = Ax[0]//100;  cB = Bx[0]//100
3235    Cen = (SGCen[cA]+SGCen[cB])%1.0
3236    cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
3237    Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
3238    Mc = np.inner(Ma,Mb.T)
3239#    print Ma,Mb,Mc
3240    Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
3241#    print Ta,Tb,Tc
3242#    print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
3243    nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
3244    #now the cell translation part
3245    if len(Ax)>1:
3246        cellA = Ax[1].split(',')
3247        cellA = [int(a) for a in cellA]
3248    else:
3249        cellA = [0,0,0]
3250    if len(Bx)>1:
3251        cellB = Bx[1].split(',')
3252        cellB = [int(b) for b in cellB]
3253    else:
3254        cellB = [0,0,0]
3255    cellC = np.add(cellA,cellB)
3256    C = str(((nC+1)+(100*cC))*(1-2*iC))+'+'+ \
3257        str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
3258    return C
3259           
3260def U2Uij(U):
3261    #returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
3262    return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
3263   
3264def Uij2U(Uij):
3265    #returns the thermal motion tensor U from Uij as numpy array
3266    return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
3267
3268def StandardizeSpcName(spcgroup):
3269    '''Accept a spacegroup name where spaces may have not been used
3270    in the names according to the GSAS convention (spaces between symmetry
3271    for each axis) and return the space group name as used in GSAS
3272    '''
3273    rspc = spcgroup.replace(' ','').upper()
3274    # deal with rhombohedral and hexagonal setting designations
3275    rhomb = ''
3276    if rspc[-1:] == 'R':
3277        rspc = rspc[:-1]
3278        rhomb = ' R'
3279    gray = ''
3280    if "1'" in rspc:
3281        gray = " 1'"
3282        rspc = rspc.replace("1'",'')
3283    elif rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
3284        rspc = rspc[:-1]
3285    else:
3286        rspc = rspc.replace("'",'')
3287    # look for a match in the spacegroup lists
3288    for i in spglist.values():
3289        for spc in i:
3290            if rspc == spc.replace(' ','').upper():
3291                return spc+gray+rhomb
3292    # how about the post-2002 orthorhombic names?
3293    if rspc in sgequiv_2002_orthorhombic:
3294        return sgequiv_2002_orthorhombic[rspc]+gray
3295    else:
3296    # not found
3297        return ''
3298
3299spgbyNum = []
3300'''Space groups indexed by number'''
3301spgbyNum = [None,
3302        'P 1','P -1',                                                   #1-2
3303        'P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
3304        'C 2/m','P 2/c','P 21/c','C 2/c',                               #3-15
3305        'P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21',
3306        'C 2 2 21','C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
3307        'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
3308        'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
3309        'C m m 2','C m c 21','C c c 2',
3310        'A m m 2','A b m 2','A m a 2','A b a 2',
3311        'F m m 2','F d d 2','I m m 2','I b a 2','I m a 2',
3312        'P m m m','P n n n','P c c m','P b a n',
3313        'P m m a','P n n a','P m n a','P c c a','P b a m','P c c n',
3314        'P b c m','P n n m','P m m n','P b c n','P b c a','P n m a',
3315        'C m c m','C m c a','C m m m','C c c m','C m m a','C c c a',
3316        'F m m m', 'F d d d',
3317        'I m m m','I b a m','I b c a','I m m a',                        #16-74
3318        'P 4','P 41','P 42','P 43',
3319        'I 4','I 41',
3320        'P -4','I -4','P 4/m','P 42/m','P 4/n','P 42/n',
3321        'I 4/m','I 41/a',
3322        'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3323        'P 42 21 2','P 43 2 2','P 43 21 2',
3324        'I 4 2 2','I 41 2 2',
3325        'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c','P 4 n c',
3326        'P 42 m c','P 42 b c',
3327        'I 4 m m','I 4 c m','I 41 m d','I 41 c d',
3328        'P -4 2 m','P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2',
3329        'P -4 c 2','P -4 b 2','P -4 n 2',
3330        'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d',
3331        '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',
3332        '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',
3333        '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',
3334        'P 42/n c m',
3335        'I 4/m m m','I 4/m c m','I 41/a m d','I 41/a c d',
3336        'P 3','P 31','P 32','R 3','P -3','R -3',
3337        'P 3 1 2','P 3 2 1','P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1',
3338        'R 3 2',
3339        'P 3 m 1','P 3 1 m','P 3 c 1','P 3 1 c',
3340        'R 3 m','R 3 c',
3341        'P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1',
3342        'R -3 m','R -3 c',                                               #75-167
3343        'P 6','P 61',
3344        'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3345        'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3346        'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3347        '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
3348        'P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m 3','P n 3',
3349        'F m -3','F d -3','I m -3',
3350        'P a -3','I a -3','P 4 3 2','P 42 3 2','F 4 3 2','F 41 3 2',
3351        'I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2','P -4 3 m',
3352        'F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
3353        'P m -3 m','P n -3 n','P m -3 n','P n -3 m',
3354        'F m -3 m','F m -3 c','F d -3 m','F d -3 c',
3355        'I m -3 m','I a -3 d',]                                       #195-230
3356spg2origins = {}
3357''' A dictionary of all spacegroups that have 2nd settings; the value is the
33581st --> 2nd setting transformation vector as X(2nd) = X(1st)-V, nonstandard ones are included.
3359'''
3360spg2origins = {
3361        'P n n n':[-.25,-.25,-.25],
3362        'P b a n':[-.25,-.25,0],'P n c b':[0,-.25,-.25],'P c n a':[-.25,0,-.25],
3363        'P m m n':[-.25,-.25,0],'P n m m':[0,-.25,-.25],'P m n m':[-.25,0,-.25],
3364        'C c c a':[0,-.25,-.25],'C c c b':[-.25,0,-.25],'A b a a':[-.25,0,-.25],
3365        'A c a a':[-.25,-.25,0],'B b c b':[-.25,-.25,0],'B b a b':[0,-.25,-.25],
3366        'F d d d':[-.125,-.125,-.125],
3367        'P 4/n':[-.25,-.25,0],'P 42/n':[-.25,-.25,-.25],'I 41/a':[0,-.25,-.125],
3368        '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],
3369        '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],
3370        'I 41/a m d':[0,.25,-.125],'I 41/a c d':[0,.25,-.125],
3371        'p n -3':[-.25,-.25,-.25],'F d -3':[-.125,-.125,-.125],'P n -3 n':[-.25,-.25,-.25],
3372        'P n -3 m':[-.25,-.25,-.25],'F d -3 m':[-.125,-.125,-.125],'F d -3 c':[-.375,-.375,-.375],
3373        'p n 3':[-.25,-.25,-.25],'F d 3':[-.125,-.125,-.125],'P n 3 n':[-.25,-.25,-.25],
3374        'P n 3 m':[-.25,-.25,-.25],'F d 3 m':[-.125,-.125,-.125],'F d - c':[-.375,-.375,-.375]}
3375spglist = {}
3376'''A dictionary of space groups as ordered and named in the pre-2002 International
3377Tables Volume A, except that spaces are used following the GSAS convention to
3378separate the different crystallographic directions.
3379Note that the symmetry codes here will recognize many non-standard space group
3380symbols with different settings. They are ordered by Laue group
3381'''
3382spglist = {
3383    'P1' : ('P 1','P -1',), # 1-2
3384    'C1' : ('C 1','C -1',),
3385    'P2/m': ('P 2','P 21','P m','P a','P c','P n',
3386        '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
3387    'C2/m':('C 2','C m','C c','C n',
3388        'C 2/m','C 2/c','C 2/n',),
3389    'A2/m':('A 2','A m','A a','A n',
3390        'A 2/m','A 2/a','A 2/n',),
3391   'Pmmm':('P 2 2 2',
3392        'P 2 2 21','P 21 2 2','P 2 21 2',
3393        'P 21 21 2','P 2 21 21','P 21 2 21',
3394        'P 21 21 21',
3395        'P m m 2','P 2 m m','P m 2 m',
3396        'P m c 21','P 21 m a','P b 21 m','P m 21 b','P c m 21','P 21 a m',
3397        'P c c 2','P 2 a a','P b 2 b',
3398        'P m a 2','P 2 m b','P c 2 m','P m 2 a','P b m 2','P 2 c m',
3399        'P c a 21','P 21 a b','P c 21 b','P b 21 a','P b c 21','P 21 c a',
3400        'P n c 2','P 2 n a','P b 2 n','P n 2 b','P c n 2','P 2 a n',
3401        'P m n 21','P 21 m n','P n 21 m','P m 21 n','P n m 21','P 21 n m',
3402        'P b a 2','P 2 c b','P c 2 a',
3403        'P n a 21','P 21 n b','P c 21 n','P n 21 a','P b n 21','P 21 c n',
3404        'P n n 2','P 2 n n','P n 2 n',
3405        'P m m m','P n n n',
3406        'P c c m','P m a a','P b m b',
3407        'P b a n','P n c b','P c n a',
3408        'P m m a','P b m m','P m c m','P m a m','P m m b','P c m m',
3409        'P n n a','P b n n','P n c n','P n a n','P n n b','P c n n',
3410        'P m n a','P b m n','P n c m','P m a n','P n m b','P c n m',
3411        'P c c a','P b a a','P b c b','P b a b','P c c b','P c a a',
3412        'P b a m','P m c b','P c m a',
3413        'P c c n','P n a a','P b n b',
3414        'P b c m','P m c a','P b m a','P c m b','P c a m','P m a b',
3415        'P n n m','P m n n','P n m n',
3416        'P m m n','P n m m','P m n m',
3417        'P b c n','P n c a','P b n a','P c n b','P c a n','P n a b',
3418        'P b c a','P c a b',
3419        'P n m a','P b n m','P m c n','P n a m','P m n b','P c m n',
3420        ),
3421    'Cmmm':('C 2 2 21','C 2 2 2','C m m 2',
3422        'C m c 21','C c m 21','C c c 2','C m 2 m','C 2 m m',
3423        'C m 2 a','C 2 m b','C c 2 m','C 2 c m','C c 2 a','C 2 c b',
3424        'C m c m','C c m m','C m c a','C c m b',
3425        'C m m m','C c c m','C m m a','C m m b','C c c a','C c c b',),
3426    'Ammm':('A 21 2 2','A 2 2 2','A 2 m m',
3427        'A 21 m a','A 21 a m','A 2 a a','A m 2 m','A m m 2',
3428        'A b m 2','A c 2 m','A m a 2','A m 2 a','A b a 2','A c 2 a',
3429        'A m m a','A m a m','A b m a','A c a m',
3430        'A m m m','A m a a','A b m m','A c m m','A c a a','A b a a',),
3431    'Bmmm':('B 2 21 2','B 2 2 2','B m 2 m',
3432        'B m 21 b','B b 21 m','B b 2 b','B m m 2','B 2 m m',
3433        'B 2 c m','B m a 2','B 2 m b','B b m 2','B 2 c b','B b a 2',
3434        'B b m m','B m m b','B b c m','B m a b',
3435        'B m m m','B b m b','B m a m','B m c m','B b a b','B b c b',),
3436    'Immm':('I 2 2 2','I 21 21 21',
3437        'I m m 2','I m 2 m','I 2 m m',
3438        'I b a 2','I 2 c b','I c 2 a',
3439        'I m a 2','I 2 m b','I c 2 m','I m 2 a','I b m 2','I 2 c m',
3440        'I m m m','I b a m','I m c b','I c m a',
3441        'I b c a','I c a b',
3442        'I m m a','I b m m ','I m c m','I m a m','I m m b','I c m m',),
3443    'Fmmm':('F 2 2 2','F m m m', 'F d d d',
3444        'F m m 2','F m 2 m','F 2 m m',
3445        'F d d 2','F d 2 d','F 2 d d',),
3446    'P4/mmm':('P 4','P 41','P 42','P 43','P -4','P 4/m','P 42/m','P 4/n','P 42/n',
3447        'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3448        'P 42 21 2','P 43 2 2','P 43 21 2','P 4 m m','P 4 b m','P 42 c m',
3449        'P 42 n m','P 4 c c','P 4 n c','P 42 m c','P 42 b c','P -4 2 m',
3450        'P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2',
3451        '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',
3452        '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',
3453        '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',
3454        'P 42/n c m',),
3455    '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',
3456        'I 4 c m','I 41 m d','I 41 c d',
3457        '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',
3458        'I 41/a m d','I 41/a c d'),
3459    'R3-H':('R 3','R -3','R 3 2','R 3 m','R 3 c','R -3 m','R -3 c',),
3460    'P6/mmm': ('P 3','P 31','P 32','P -3','P 3 1 2','P 3 2 1','P 31 1 2',
3461        'P 31 2 1','P 32 1 2','P 32 2 1', 'P 3 m 1','P 3 1 m','P 3 c 1',
3462        'P 3 1 c','P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1','P 6','P 61',
3463        'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3464        'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3465        'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3466        'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',),
3467    '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',
3468        'P 4 3 2','P 42 3 2','P 43 3 2','P 41 3 2','P -4 3 m','P -4 3 n',
3469        'P m 3 m','P m -3 m','P n 3 n','P n -3 n',
3470        'P m 3 n','P m -3 n','P n 3 m','P n -3 m',),
3471    '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',
3472        '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'),
3473    'Fm3m':('F 2 3','F m 3','F m -3','F d 3','F d -3',
3474        'F 4 3 2','F 41 3 2','F -4 3 m','F -4 3 c',
3475        'F m 3 m','F m -3 m','F m 3 c','F m -3 c',
3476        'F d 3 m','F d -3 m','F d 3 c','F d -3 c',),
3477}
3478sgequiv_2002_orthorhombic = {}
3479''' A dictionary of orthorhombic space groups that were renamed in the 2002 Volume A,
3480 along with the pre-2002 name. The e designates a double glide-plane
3481'''
3482sgequiv_2002_orthorhombic = {
3483        'AEM2':'A b m 2','B2EM':'B 2 c m','CM2E':'C m 2 a',
3484        'AE2M':'A c 2 m','BME2':'B m a 2','C2ME':'C 2 m b',
3485        'AEA2':'A b a 2','B2EB':'B 2 c b','CC2E':'C c 2 a',
3486        'AE2A':'A c 2 a','BBE2':'B b a 2','C2CE':'C 2 c b',
3487        'CMCE':'C m c a','AEMA':'A b m a','BBEM':'B b c m',
3488        'BMEB':'B m a b','CCME':'C c m b','AEAM':'A c a m',
3489        'CMME':'C m m a','AEMM':'A b m m','BMEM':'B m c m',
3490        'CCCE':'C c c a','AEAA':'A b a a','BBEB':'B b c b'}
3491
3492#'A few non-standard space groups for test use'
3493nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h', 
3494                      'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
3495                      'R 3 c r','R -3 c r','R -3 m r',),
3496
3497#Use the space groups types in this order to list the symbols in the
3498#order they are listed in the International Tables, vol. A'''
3499symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 
3500               'trigonal', 'hexagonal', 'cubic')
3501
3502# self-test materials follow. Requires files in directory testinp
3503selftestlist = []
3504'''Defines a list of self-tests'''
3505selftestquiet = True
3506def _ReportTest():
3507    'Report name and doc string of current routine when ``selftestquiet`` is False'
3508    if not selftestquiet:
3509        import inspect
3510        caller = inspect.stack()[1][3]
3511        doc = eval(caller).__doc__
3512        if doc is not None:
3513            print('testing '+__file__+' with '+caller+' ('+doc+')')
3514        else:
3515            print('testing '+__file__()+" with "+caller)
3516def test0():
3517    '''self-test #0: exercise MoveToUnitCell'''
3518    _ReportTest()
3519    msg = "MoveToUnitCell failed"
3520    assert (MoveToUnitCell([1,2,3])[0] == [0,0,0]).all, msg
3521    assert (MoveToUnitCell([2,-1,-2])[0] == [0,0,0]).all, msg
3522    assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9)[0] < 1e-6, msg
3523    assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1)[0] < 1e-6, msg
3524selftestlist.append(test0)
3525
3526def test1():
3527    '''self-test #1: SpcGroup against previous results'''
3528    #'''self-test #1: SpcGroup and SGPrint against previous results'''
3529    _ReportTest()
3530    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3531    if ospath.exists(testdir):
3532        if testdir not in sys.path: sys.path.insert(0,testdir)
3533    import spctestinp
3534    def CompareSpcGroup(spc, referr, refdict, reflist): 
3535        'Compare output from GSASIIspc.SpcGroup with results from a previous run'
3536        # if an error is reported, the dictionary can be ignored
3537        msg0 = "CompareSpcGroup failed on space group %s" % spc
3538        result = SpcGroup(spc)
3539        if result[0] == referr and referr > 0: return True
3540#        #print result[1]['SpGrp']
3541        #msg = msg0 + " in list lengths"
3542        #assert len(keys) == len(refdict.keys()), msg
3543