source: trunk/GSASIIspc.py @ 3379

Last change on this file since 3379 was 3379, checked in by vondreele, 5 years ago

fix problem of bad space group symbol (no field spaces) found in cif file

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