source: trunk/GSASIIspc.py @ 3307

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

provide print tool for magnetic space group operators that can be fed into Bilbao routine for determinimg std BNS space group with transformation

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