source: trunk/GSASIIspc.py @ 3320

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

implement new delt/sig subplot for PWDR & SASD plots - nvoked with 'w' key
make new MagStructureFactor2 routine - removed magnetism stuff from StructureFactor2

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