source: trunk/GSASIIspc.py @ 3317

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

magnetic site symmetry restrictions on moments & magnetic site symmetry symbols

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