source: trunk/GSASIIspc.py @ 3308

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

more fixes for mag site symmetry constraints

  • 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-05 21:31:19 +0000 (Mon, 05 Mar 2018) $
12# $Author: vondreele $
13# $Revision: 3308 $
14# $URL: trunk/GSASIIspc.py $
15# $Id: GSASIIspc.py 3308 2018-03-05 21:31:19Z 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: 3308 $")
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']:
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','m3','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        SGData['SGSpin'][-1] = -1
793        if 'R' in BNS:
794            SGData['SGSpin'][-1] = -1
795    elif '(S)' in BNS:
796        SGData['SGSpin'][-1] = -1
797        SGData['SGSpin'] += [-1,-1,-1,]
798        Tmat *= 2.0
799    SGData['SGSpin'].append(-1)
800    if 'P' not in BNS:
801        SGData['SGSpin'].append(-1)
802    C = SGCen+A
803    SGData['SGCen'] = np.vstack((SGCen,C))%1.
804    return Tmat
805   
806def CheckSpin(isym,SGData):
807    ''' Check for exceptions in spin rules
808    '''
809    if SGData['SpGrp'] in ['C c','C 1 c 1','A a','A 1 a 1','B b 1 1','C c 1 1',
810        'A 1 1 a','B 1 1 b','I -4']:
811        if SGData['SGSpin'][:2] == [-1,-1]:
812            SGData['SGSpin'][(isym+1)%2] = 1
813    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',
814        'A 1 1 2/a','B 1 1 2/b']:
815        if SGData['SGSpin'][1:3] == [-1,-1]:
816            SGData['SGSpin'][isym%2+1] = 1
817    elif SGData['SGPtGrp'] in ['222','mm2','2mm','m2m']:
818        if SGData['SGSpin'][0]*SGData['SGSpin'][1]*SGData['SGSpin'][2] < 0:
819            SGData['SGSpin'][(isym+1)%3] *= -1
820        if SGData['SpGrp'][0] == 'F' and isym > 2:
821            SGData['SGSpin'][(isym+1)%3+3] *= -1
822    elif SGData['SGPtGrp'] == 'mmm':
823        if SGData['SpGrp'][0] == 'F' and isym > 2:
824            SGData['SGSpin'][(isym+1)%3+3] *= -1
825        elif SGData['SGSpin'][3] < 0:
826            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',
827                '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']:
828                for i in [0,1,2]:
829                    if i != isym and SGData['SGSpin'][i] < 0:
830                        SGData['SGSpin'][i] = 1
831            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']:
832                if SGData['SGSpin'][0]*SGData['SGSpin'][1]*SGData['SGSpin'][2] < 0:
833                    SGData['SGSpin'][(isym+1)%3] *= -1
834    elif SGData['SpGrp'] in ['I -4 m 2','I -4 c 2']:
835        if SGData['SGSpin'][2] < 0:
836            if 'm' in SGData['SpGrp']:
837                SGData['SGSpin'][1] = 1
838            elif isym < 2:
839                if SGData['SGSpin'][isym] < 0:
840                    SGData['SGSpin'][:2] = [-1,-1]
841                else:
842                    SGData['SGSpin'][:2] = [1,1]
843            else:
844                SGData['SGSpin'][:2] = [1,1]
845   
846def MagSGSym(SGData):       #needs to use SGPtGrp not SGLaue!
847    SGLaue = SGData['SGLaue']
848    if '1' not in SGData['GenSym']:        #patch for old gpx files
849        SGData['GenSym'] = ['1',]+SGData['GenSym']
850        SGData['SGSpin'] = [1,]+list(SGData['SGSpin'])
851    if len(SGData['SGSpin'])<len(SGData['GenSym']):
852        SGData['SGSpin'] = [1,]+list(SGData['SGSpin'])      #end patch
853    GenSym = SGData['GenSym'][1:]       #skip identity
854    SpnFlp = SGData['SGSpin']
855#    print('SpnFlp',SpnFlp)
856    SGPtGrp = SGData['SGPtGrp']
857    if len(SpnFlp) == 1:
858        SGData['MagPtGp'] = SGPtGrp
859        return SGData['SpGrp']
860    magSym = SGData['SpGrp'].split()
861    if SGLaue in ['-1',]:
862        SGData['MagPtGp'] = SGPtGrp
863        if SpnFlp[1] == -1:
864            magSym[1] += "'"
865            SGData['MagPtGp'] += "'"
866    elif SGLaue in ['2/m','4/m','6/m']: #all ok
867        Uniq = {'a':1,'b':2,'c':3,'':1}
868        id = [0,1]
869        if len(magSym) > 2:
870            id = [0,Uniq[SGData['SGUniq']]]
871        sym = magSym[id[1]].split('/')
872        Ptsym = SGLaue.split('/')
873        if len(GenSym) == 3:
874            for i in [0,1,2]:
875                if SpnFlp[i+1] < 0:
876                    sym[i] += "'"
877                    Ptsym[i] += "'"
878        else:
879            for i in range(len(GenSym)):
880                if SpnFlp[i+1] < 0:                     
881                    sym[i] += "'"
882                    Ptsym[i] += "'"
883        SGData['MagPtGp'] = '/'.join(Ptsym)
884        magSym[id[1]] = '/'.join(sym)
885    elif SGPtGrp in ['mmm','mm2','m2m','2mm','222']:
886        SGData['MagPtGp'] = ''
887        for i in [0,1,2]:
888            SGData['MagPtGp'] += SGPtGrp[i]
889            if SpnFlp[i+1] < 0:
890                magSym[i+1] += "'"
891                SGData['MagPtGp'] += "'"
892    elif SGLaue == '6/mmm': #ok
893        magPtGp = list(SGPtGrp)
894        if len(GenSym) == 2:
895            for i in [0,1]:
896                if SpnFlp[i+1] < 0:
897                    magSym[i+2] += "'"
898                    magPtGp[i+1] += "'"
899            if SpnFlp[1]*SpnFlp[2] < 0:
900                magSym[1] += "'"
901                magPtGp[0] += "'"
902        else:
903            sym = magSym[1].split('/')
904            Ptsym = ['6','m']
905            magPtGp = ['','m','m']
906            for i in [0,1,2]:
907                if SpnFlp[i+1] < 0:
908                    if i:
909                        magSym[i+1] += "'"
910                        magPtGp[i] += "'"
911                    else:
912                        sym[1] += "'"
913                        Ptsym[0] += "'"
914            if SpnFlp[2]*SpnFlp[3] < 0:
915                sym[0] += "'"                   
916                Ptsym[0] += "'"                   
917            magSym[1] = '/'.join(sym)
918            magPtGp[0] = '/'.join(Ptsym)
919        SGData['MagPtGp'] = ''.join(magPtGp)
920    elif SGLaue == '4/mmm':
921        magPtGp = list(SGPtGrp)
922        if len(GenSym) == 2:
923            for i in [0,1]:
924                if SpnFlp[i+1] < 0:
925                    magSym[i+2] += "'"
926                    magPtGp[i+1] += "'"
927            if SpnFlp[1]*SpnFlp[2] < 0:
928                magSym[1] += "'"
929                magPtGp[0] += "'"
930        else:
931            if '/' in magSym[1]:    #P 4/m m m, etc.
932                sym = magSym[1].split('/')
933                Ptsym = ['4','m']
934                magPtGp = ['','m','m']
935                for i in [0,1,2]:
936                    if SpnFlp[i+1] < 0:
937                        if i:
938                            magSym[i+1] += "'"
939                            magPtGp[i] += "'"
940                        else:
941                            sym[1] += "'"
942                            Ptsym[1] += "'"
943                if SpnFlp[2]*SpnFlp[3] < 0:
944                    sym[0] += "'"                   
945                    Ptsym[0] += "'"                   
946                magSym[1] = '/'.join(sym)
947                magPtGp[0] = '/'.join(Ptsym)
948            else:
949                for i in [0,1]:
950                    if SpnFlp[i+1] < 0:
951                        magSym[i+2] += "'"
952                if SpnFlp[1]*SpnFlp[2] < 0:
953                    magSym[1] += "'"
954        SGData['MagPtGp'] = ''.join(magPtGp)
955    elif SGLaue in ['3','3m1','31m']:   #ok
956        Ptsym = list(SGLaue)
957        if len(GenSym) == 1:    #all ok
958            id = 2
959            if (len(magSym) == 4) and (magSym[2] == '1'):
960                id = 3
961            if '3' in GenSym[0]:
962                id = 1
963            magSym[id].strip("'")
964            if SpnFlp[1] < 0:
965                magSym[id] += "'"
966                Ptsym[id-1] += "'"
967        elif len(GenSym) == 2:
968            if 'R' in GenSym[1]:
969                magSym[-1].strip("'")
970                if SpnFlp[1] < 0:
971                    magSym[-1] += "'"
972                    Ptsym[-1] += "'"
973            else:
974                i,j = [1,2]
975                if magSym[2] == '1':
976                    i,j = [1,3]
977                magSym[i].strip("'")
978                Ptsym[i-1].strip("'")
979                magSym[j].strip("'")
980                Ptsym[j-1].strip("'")
981                if SpnFlp[1:3] == [1,-1]:
982                    magSym[i] += "'"
983                    Ptsym[i-1] += "'"
984                elif SpnFlp[1:3] == [-1,-1]:
985                    magSym[j] += "'"
986                    Ptsym[j-1] += "'"
987                elif SpnFlp[1:3] == [-1,1]:
988                    magSym[i] += "'"
989                    Ptsym[i-1] += "'"
990                    magSym[j] += "'"
991                    Ptsym[j-1] += "'"
992        elif len(GenSym):
993            if 'c' not in magSym[2]:
994                i,j = [1,2]
995                magSym[i].strip("'")
996                Ptsym[i-1].strip("'")
997                magSym[j].strip("'")
998                Ptsym[j-1].strip("'")
999                if SpnFlp[1:3] == [1,-1]:
1000                    magSym[i] += "'"
1001                    Ptsym[i-1] += "'"
1002                elif SpnFlp[1:3] == [-1,-1]:
1003                    magSym[j] += "'"
1004                    Ptsym[j-1] += "'"
1005                elif SpnFlp[2] == [-1,1]:
1006                    magSym[i] += "'"
1007                    Ptsym[i-1] += "'"
1008                    magSym[j] += "'"
1009                    Ptsym[j-1] += "'"
1010        SGData['MagPtGp'] = ''.join(Ptsym)
1011    elif SGData['SGPtGrp'] == '23' and len(magSym):
1012        SGData['MagPtGp'] = '23'
1013    elif SGData['SGPtGrp'] == 'm3':
1014        SGData['MagPtGp'] = "m3"
1015        if SpnFlp[1] < 0:
1016            magSym[1] += "'"
1017            magSym[2] += "'"
1018            SGData['MagPtGp'] = "m'3'"
1019        if SpnFlp[2] < 0:
1020            if not 'm' in magSym[1]:    #only Ia3
1021                magSym[1].strip("'")
1022                SGData['MagPtGp'] = "m3'"
1023    elif SGData['SGPtGrp'] in ['432','-43m']:
1024        Ptsym = SGData['SGPtGrp'].split('3')
1025        if SpnFlp[1] < 0:
1026            magSym[1] += "'"
1027            Ptsym[0] += "'"
1028            magSym[3] += "'"
1029            Ptsym[1] += "'"
1030        SGData['MagPtGp'] = '3'.join(Ptsym)
1031    elif SGData['SGPtGrp'] == 'm-3m':
1032        Ptsym = ['m','3','m']
1033        if SpnFlp[1:3] == [-1,1]:
1034            magSym[1] += "'"
1035            Ptsym[0] += "'"
1036            magSym[2] += "'"
1037            Ptsym[1] += "'"
1038        elif SpnFlp[1:3] == [1,-1]:
1039            magSym[3] += "'"
1040            Ptsym[2] += "'"
1041        elif SpnFlp[1:3] == [-1,-1]:
1042            magSym[1] += "'"
1043            Ptsym[0] += "'"
1044            magSym[2] += "'"
1045            Ptsym[1] += "'"
1046            magSym[3] += "'"
1047            Ptsym[2] += "'"
1048        SGData['MagPtGp'] = ''.join(Ptsym)
1049#    print SpnFlp
1050    magSym[0] = SGData.get('BNSlattsym',[SGData['SGLatt'],[0,0,0]])[0]
1051    return ' '.join(magSym)
1052   
1053def MagText2MTS(mcifOpr,CIF=True):
1054    "From magnetic space group cif text returns matrix/translation + spin flip"
1055    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],
1056           'z':[0,0,1],'+z':[0,0,1],'-z':[0,0,-1],'x-y':[1,-1,0],'-x+y':[-1,1,0],'y-x':[-1,1,0]}
1057    ops = mcifOpr.split(",")
1058    M = []
1059    T = []
1060    for op in ops[:3]:
1061        ip = len(op)
1062        if '/' in op:
1063            if CIF:
1064                opMT = op.split('+')
1065                T.append(eval(opMT[1]))
1066            else:
1067                ip = op.index('/')
1068                T.append(eval(op[:ip+2]))
1069                opMT = [op[ip+2:],'']
1070        else:
1071            opMT = [op,'']
1072            T.append(0.)
1073        M.append(XYZ[opMT[0].lower()])
1074    spnflp = 1
1075    if '-1' in ops[3]:
1076        spnflp = -1
1077    return np.array(M),np.array(T),spnflp
1078           
1079def MagSSText2MTS(mcifOpr):
1080    "From magnetic super space group cif text returns matrix/translation + spin flip"
1081    XYZ = {'x1':[1,0,0,0],'-x1':[-1,0,0,0],
1082           'x2':[0,1,0,0],'-x2':[0,-1,0,0],
1083           'x3':[0,0,1,0],'-x3':[0,0,-1,0],
1084           'x4':[0,0,0,1],'-x4':[0,0,0,-1],
1085           'x1-x2':[1,-1,0,0],'-x1+x2':[-1,1,0,0],
1086           'x1-x4':[1,0,0,-1],'-x1+x4':[-1,0,0,1],
1087           'x2-x4':[0,1,0,-1],'-x2+x4':[0,-1,0,1],
1088           '-x1-x2+x4':[-1,-1,0,1],'x1+x2-x4':[1,1,0,-1]}
1089    ops = mcifOpr.split(",")
1090    M = []
1091    T = []
1092    for op in ops[:4]:
1093        ip = len(op)
1094        if '/' in op:
1095            ip = op.index('/')-2
1096            T.append(eval(op[ip:]))
1097        else:
1098            T.append(0.)
1099        M.append(XYZ[op[:ip]])
1100    spnflp = 1
1101    if '-1' in ops[4]:
1102        spnflp = -1
1103    return np.array(M),np.array(T),spnflp
1104           
1105def GenMagOps(SGData):
1106    FlpSpn = SGData['SGSpin']
1107    Nsym = len(SGData['SGOps'])
1108    Ncv = len(SGData['SGCen'])
1109    sgOp = [M for M,T in SGData['SGOps']]
1110    OprName = [GetOprPtrName(str(irtx))[1] for irtx in PackRot(SGData['SGOps'])]
1111    if SGData['SGInv'] and not SGData.get('SGFixed',False):
1112        Nsym *= 2
1113        sgOp += [-M for M,T in SGData['SGOps']]
1114        OprName += [GetOprPtrName(str(-irtx))[1] for irtx in PackRot(SGData['SGOps'])]
1115    Nsyms = 0
1116    sgOps = []
1117    OprNames = []
1118    for incv in range(Ncv):
1119        Nsyms += Nsym
1120        sgOps += sgOp
1121        OprNames += OprName   
1122    SpnFlp = np.ones(Nsym,dtype=np.int)
1123    GenFlg = SGData.get('GenFlg',[0])
1124#    print ('GenFlg:',SGData['GenFlg'])
1125#    print ('GenSym:',SGData['GenSym'])
1126    Nfl = len(GenFlg)
1127    for ieqv in range(Nsym):
1128        for iunq in range(Nfl):
1129            if SGData['SGGen'][ieqv] & GenFlg[iunq]:
1130                SpnFlp[ieqv] *= FlpSpn[iunq]
1131#        print ('\nMagSpGrp:',SGData['MagSpGrp'],Ncv)
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    print('for ',siteSym)
2193    for opr in dupDir:
2194        indx = GetNXUPQsym(opr)
2195        if SpnFlp[dupDir[opr]] > 0.:
2196            csi = CSxinel[indx[2]]  #P
2197        else:
2198            csi = CSxinel[indx[3]]  #Q
2199        print(opr,SpnFlp[dupDir[opr]],indx,csi,CSI)
2200        if not len(csi):
2201            return [[0,0,0],[0.,0.,0.]]
2202        for kcs in [0,1,2]:
2203            if csi[0][kcs] == 0 and CSI[0][kcs] != 0:
2204                jcs = CSI[0][kcs]
2205                for ics in [0,1,2]:
2206                    if CSI[0][ics] == jcs:
2207                        CSI[0][ics] = 0
2208                        CSI[1][ics] = 0.
2209                    elif CSI[0][ics] > jcs:
2210                        CSI[0][ics] = CSI[0][ics]-1
2211            elif (CSI[0][kcs] == csi[0][kcs]) and (CSI[1][kcs] != csi[1][kcs]):
2212                CSI[1][kcs] = csi[1][kcs]
2213            elif CSI[0][kcs] >= csi[0][kcs]:
2214                CSI[0][kcs] = min(CSI[0][kcs],csi[0][kcs])
2215                if CSI[1][kcs] != csi[1][kcs]:
2216                    if CSI[1][kcs] == 1.:
2217                        CSI[1][kcs] = csi[1][kcs]
2218        print(CSI)
2219    return CSI
2220   
2221def getTauT(tau,sop,ssop,XYZ,wave=np.zeros(3)):
2222    phase = np.sum(XYZ*wave)
2223    ssopinv = nl.inv(ssop[0])
2224    mst = ssopinv[3][:3]
2225    epsinv = ssopinv[3][3]
2226    sdet = nl.det(sop[0])
2227    ssdet = nl.det(ssop[0])
2228    dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
2229    dT = 1.0
2230    if np.any(dtau%.5):
2231        sumdtau = np.sum(dtau%.5)
2232        dT = 0.
2233        if np.abs(sumdtau-.5) > 1.e-4:
2234            dT = np.tan(np.pi*sumdtau)
2235    tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]+phase)
2236    return sdet,ssdet,dtau,dT,tauT
2237   
2238def OpsfromStringOps(A,SGData,SSGData):
2239    SGOps = SGData['SGOps']
2240    SSGOps = SSGData['SSGOps']
2241    Ax = A.split('+')
2242    Ax[0] = int(Ax[0])
2243    iC = 1
2244    if Ax[0] < 0:
2245        iC = -1
2246    Ax[0] = abs(Ax[0])
2247    nA = Ax[0]%100-1
2248    nC = Ax[0]//100
2249    unit = [0,0,0]
2250    if len(Ax) > 1:
2251        unit = eval('['+Ax[1]+']')
2252    return SGOps[nA],SSGOps[nA],iC,SGData['SGCen'][nC],unit
2253   
2254def GetSSfxuinel(waveType,Stype,nH,XYZ,SGData,SSGData,debug=False):
2255   
2256    def orderParms(CSI):
2257        parms = [0,]
2258        for csi in CSI:
2259            for i in [0,1,2]:
2260                if csi[i] not in parms:
2261                    parms.append(csi[i])
2262        for csi in CSI:
2263            for i in [0,1,2]:
2264                csi[i] = parms.index(csi[i])
2265        return CSI
2266       
2267    def fracCrenel(tau,Toff,Twid):
2268        Tau = (tau-Toff[:,np.newaxis])%1.
2269        A = np.where(Tau<Twid[:,np.newaxis],1.,0.)
2270        return A
2271       
2272    def fracFourier(tau,nH,fsin,fcos):
2273        SA = np.sin(2.*nH*np.pi*tau)
2274        CB = np.cos(2.*nH*np.pi*tau)
2275        A = SA[np.newaxis,np.newaxis,:]*fsin[:,:,np.newaxis]
2276        B = CB[np.newaxis,np.newaxis,:]*fcos[:,:,np.newaxis]
2277        return A+B
2278       
2279    def posFourier(tau,nH,psin,pcos):
2280        SA = np.sin(2*nH*np.pi*tau)
2281        CB = np.cos(2*nH*np.pi*tau)
2282        A = SA[np.newaxis,np.newaxis,:]*psin[:,:,np.newaxis]
2283        B = CB[np.newaxis,np.newaxis,:]*pcos[:,:,np.newaxis]
2284        return A+B   
2285
2286    def posSawtooth(tau,Toff,slopes):
2287        Tau = (tau-Toff)%1.
2288        A = slopes[:,np.newaxis]*Tau
2289        return A
2290   
2291    def posZigZag(tau,Tmm,XYZmax):
2292        DT = Tmm[1]-Tmm[0]
2293        slopeUp = 2.*XYZmax/DT
2294        slopeDn = 2.*XYZmax/(1.-DT)
2295        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])
2296        return A
2297
2298    def posBlock(tau,Tmm,XYZmax):
2299        A = np.array([np.where(Tmm[0] < t <= Tmm[1],XYZmax,-XYZmax) for t in tau])
2300        return A
2301       
2302    def DoFrac():
2303        dF = None
2304        dFTP = None
2305        if siteSym == '1':
2306            CSI = [[1,0],[2,0]],2*[[1.,0.],]
2307        elif siteSym == '-1':
2308            CSI = [[1,0],[0,0]],2*[[1.,0.],]
2309        else:
2310            delt2 = np.eye(2)*0.001
2311            FSC = np.ones(2,dtype='i')
2312            CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
2313            if 'Crenel' in waveType:
2314                dF = np.zeros_like(tau)
2315            else:
2316                dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2317            dFT = np.zeros_like(dF)
2318            dFTP = []
2319            for i in SdIndx:
2320                sop = Sop[i]
2321                ssop = SSop[i]           
2322                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2323                fsc = np.ones(2,dtype='i')
2324                if 'Crenel' in waveType:
2325                    dFT = np.zeros_like(tau)
2326                    fsc = [1,1]
2327                else:   #Fourier
2328                    dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
2329                    dFT = nl.det(sop[0])*dFT
2330                    dFT = dFT[:,np.argsort(tauT)]
2331                    dFT[0] *= ssdet
2332                    dFT[1] *= sdet
2333                    dFTP.append(dFT)
2334               
2335                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2336                        fsc = [1,1]
2337                        if dT:
2338                            CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
2339                        else:
2340                            CSI = [[[1,0],[0,0]],[[1.,0.],[0.,0.]]]
2341                        FSC = np.zeros(2,dtype='i')
2342                        return CSI,dF,dFTP
2343                    else:
2344                        for i in range(2):
2345                            if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
2346                                fsc[i] = 1
2347                            else:
2348                                fsc[i] = 0
2349                        FSC &= fsc
2350                        if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc)
2351            n = -1
2352            for i,F in enumerate(FSC):
2353                if F:
2354                    n += 1
2355                    CSI[0][i] = n+1
2356                    CSI[1][i] = 1.0
2357           
2358        return CSI,dF,dFTP
2359       
2360    def DoXYZ():
2361        dX,dXTP = None,None
2362        if siteSym == '1':
2363            CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2364        elif siteSym == '-1':
2365            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.],]
2366        else:
2367            delt4 = np.ones(4)*0.001
2368            delt5 = np.ones(5)*0.001
2369            delt6 = np.eye(6)*0.001
2370            if 'Fourier' in waveType:
2371                dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2372                  #3x6x12 modulated position array (X,Spos,tau)& force positive
2373                CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2374            elif waveType == 'Sawtooth':
2375                dX = posSawtooth(tau,delt4[0],delt4[1:])
2376                CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0]]),
2377                    np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2378            elif waveType in ['ZigZag','Block']:
2379                if waveType == 'ZigZag':
2380                    dX = posZigZag(tau,delt5[:2],delt5[2:])
2381                else:
2382                    dX = posBlock(tau,delt5[:2],delt5[2:])
2383                CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2384                    np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2385            XSC = np.ones(6,dtype='i')
2386            dXTP = []
2387            for i in SdIndx:
2388                sop = Sop[i]
2389                ssop = SSop[i]
2390                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2391                xsc = np.ones(6,dtype='i')
2392                if 'Fourier' in waveType:
2393                    dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2394                elif waveType == 'Sawtooth':
2395                    dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2396                elif waveType == 'ZigZag':
2397                    dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2398                elif waveType == 'Block':
2399                    dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2400                dXT = np.inner(sop[0],dXT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2401                dXT = np.swapaxes(dXT,1,2)      # back to array(3x6x49)
2402                dXT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2403                dXTP.append(dXT)
2404                if waveType == 'Fourier':
2405                    for i in range(3):
2406                        if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2407                            xsc[i] = 0
2408                        if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2409                            xsc[i+3] = 0
2410                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2411                        xsc[3:6] = 0
2412                        CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2413                            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2414                        if dT:
2415                            if '(x)' in siteSym:
2416                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2417                                if 'm' in siteSym and len(SdIndx) == 1:
2418                                    CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2419                            elif '(y)' in siteSym:
2420                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2421                                if 'm' in siteSym and len(SdIndx) == 1:
2422                                    CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2423                            elif '(z)' in siteSym:
2424                                CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2425                                if 'm' in siteSym and len(SdIndx) == 1:
2426                                    CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2427                        else:
2428                            CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2429                    if '4/mmm' in laue:
2430                        if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2431                            if '(xy)' in siteSym:
2432                                CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2433                                if dT:
2434                                    CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2435                                else:
2436                                    CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2437                        if '(xy)' in siteSym or '(+-0)' in siteSym:
2438                            mul = 1
2439                            if '(+-0)' in siteSym:
2440                                mul = -1
2441                            if np.allclose(dX[0,0,:],dXT[1,0,:]):
2442                                CSI[0][3:5] = [[11,0,0],[11,0,0]]
2443                                CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2444                                xsc[3:5] = 0
2445                            if np.allclose(dX[0,3,:],dXT[0,4,:]):
2446                                CSI[0][:2] = [[12,0,0],[12,0,0]]
2447                                CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2448                                xsc[:2] = 0
2449                XSC &= xsc
2450                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2451            if waveType == 'Fourier':
2452                n = -1
2453                if debug: print (XSC)
2454                for i,X in enumerate(XSC):
2455                    if X:
2456                        n += 1
2457                        CSI[0][i][0] = n+1
2458                        CSI[1][i][0] = 1.0
2459           
2460        return list(CSI),dX,dXTP
2461       
2462    def DoUij():
2463        dU,dUTP = None,None
2464        if siteSym == '1':
2465            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2466                [7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2467        elif siteSym == '-1':
2468            CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],   \
2469                6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2470        else:
2471            tau = np.linspace(0,1,49,True)
2472            delt12 = np.eye(12)*0.0001
2473            dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
2474            CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2475            USC = np.ones(12,dtype='i')
2476            dUTP = []
2477            dtau = 0.
2478            for i in SdIndx:
2479                sop = Sop[i]
2480                ssop = SSop[i]
2481                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2482                usc = np.ones(12,dtype='i')
2483                dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
2484                dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3)    #convert dUT to 12x49x3x3
2485                dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2486                dUT = np.array(U2Uij(dUijT))    #convert to 6x12x49
2487                dUT = dUT[:,:,np.argsort(tauT)]
2488                dUT[:,:6,:] *=(ssdet*sdet)
2489                dUTP.append(dUT)
2490                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2491                    if dT:
2492                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2493                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2494                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2495                        [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2496                    else:
2497                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2498                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2499                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2500                        [0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2501                    if 'mm2(x)' in siteSym and dT:
2502                        CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2503                        USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2504                    elif '(xy)' in siteSym and dT:
2505                        CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2506                            [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2507                        CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2508                        USC = [1,1,1,1,1,1,1,1,1,1,1,1]                             
2509                    elif '(x)' in siteSym and dT:
2510                        CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2511                    elif '(y)' in siteSym and dT:
2512                        CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2513                    elif '(z)' in siteSym and dT:
2514                        CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2515                    for i in range(6):
2516                        if not USC[i]:
2517                            CSI[0][i] = [0,0,0]
2518                            CSI[1][i] = [0.,0.,0.]
2519                            CSI[0][i+6] = [0,0,0]
2520                            CSI[1][i+6] = [0.,0.,0.]
2521                else:                       
2522                    for i in range(6):
2523                        if not np.allclose(dU[i,i,:],dUT[i,i,:]):  #sin part
2524                            usc[i] = 0
2525                        if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]):   #cos part
2526                            usc[i+6] = 0
2527                    if np.any(dUT[1,0,:]):
2528                        if '4/m' in siteSym:
2529                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2530                            if ssop[1][3]:
2531                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2532                                usc[9] = 1
2533                            else:
2534                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2535                                usc[9] = 0
2536                        elif '4' in siteSym:
2537                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2538                            CSI[0][:2] = [[11,0,0],[11,0,0]]
2539                            if ssop[1][3]:
2540                                CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2541                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2542                                usc[2] = 0
2543                                usc[8] = 0
2544                                usc[3] = 1
2545                                usc[9] = 1
2546                            else:
2547                                CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2548                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2549                                usc[2] = 1
2550                                usc[8] = 1
2551                                usc[3] = 0               
2552                                usc[9] = 0
2553                        elif 'xy' in siteSym or '+-0' in siteSym:
2554                            if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2555                                CSI[0][4:6] = [[12,0,0],[12,0,0]]
2556                                CSI[0][6:8] = [[11,0,0],[11,0,0]]
2557                                CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2558                                CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2559                                usc[4:6] = 0
2560                                usc[6:8] = 0
2561                           
2562                    if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2563                USC &= usc
2564            if debug: print (USC)
2565            if not np.any(dtau%.5):
2566                n = -1
2567                for i,U in enumerate(USC):
2568                    if U:
2569                        n += 1
2570                        CSI[0][i][0] = n+1
2571                        CSI[1][i][0] = 1.0
2572   
2573        return list(CSI),dU,dUTP
2574   
2575    def DoMag():
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        if siteSym == '1':
2578            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2579        elif siteSym in ['-1','mmm','2/m(x)','2/m(y)','2/m(z)','4/mmm001']:
2580            CSI = 3*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3*[[0.,0.,0.],]+3*[[1.,0.,0.],]
2581        else:
2582            delt6 = np.eye(6)*0.001
2583            dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,np.newaxis,np.newaxis]
2584              #3x6x12 modulated moment array (M,Spos,tau)& force positive
2585            CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2586            MSC = np.ones(6,dtype='i')
2587            dMTP = []
2588            for i in SdIndx:
2589                sop = Sop[i]
2590                ssop = SSop[i]
2591                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2592                msc = np.ones(6,dtype='i')
2593                dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2594                dMT = np.inner(sop[0],dMT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2595                dMT = np.swapaxes(dMT,1,2)      # back to array(3x6x49)
2596                dMT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2597                dMTP.append(dMT)
2598                for i in range(3):
2599                    if not np.allclose(dM[i,i,:],dMT[i,i,:]):
2600                        msc[i] = 0
2601                    if not np.allclose(dM[i,i+3,:],dMT[i,i+3,:]):
2602                        msc[i+3] = 0
2603                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2604                    msc[3:6] = 0
2605                    CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2606                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2607                    if dT:
2608                        if '(x)' in siteSym:
2609                            CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2610                            if 'm' in siteSym and len(SdIndx) == 1:
2611                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2612                        elif '(y)' in siteSym:
2613                            CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2614                            if 'm' in siteSym and len(SdIndx) == 1:
2615                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2616                        elif '(z)' in siteSym:
2617                            CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2618                            if 'm' in siteSym and len(SdIndx) == 1:
2619                                CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2620                    else:
2621                        CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2622                if '4/mmm' in laue:
2623                    if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2624                        if '(xy)' in siteSym:
2625                            CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2626                            if dT:
2627                                CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2628                            else:
2629                                CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2630                    if '(xy)' in siteSym or '(+-0)' in siteSym:
2631                        mul = 1
2632                        if '(+-0)' in siteSym:
2633                            mul = -1
2634                        if np.allclose(dM[0,0,:],dMT[1,0,:]):
2635                            CSI[0][3:5] = [[11,0,0],[11,0,0]]
2636                            CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2637                            msc[3:5] = 0
2638                        if np.allclose(dM[0,3,:],dMT[0,4,:]):
2639                            CSI[0][:2] = [[12,0,0],[12,0,0]]
2640                            CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2641                            msc[:2] = 0
2642                MSC &= msc
2643                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2644            n = -1
2645            if debug: print (MSC)
2646            for i,M in enumerate(MSC):
2647                if M:
2648                    n += 1
2649                    CSI[0][i][0] = n+1
2650                    CSI[1][i][0] = 1.0
2651
2652        return CSI,None,None
2653       
2654    if debug: print ('super space group: '+SSGData['SSpGrp'])
2655    xyz = np.array(XYZ)%1.
2656    SGOps = copy.deepcopy(SGData['SGOps'])
2657    laue = SGData['SGLaue']
2658    siteSym = SytSym(XYZ,SGData)[0].strip()
2659    if debug: print ('siteSym: '+siteSym)
2660    SSGOps = copy.deepcopy(SSGData['SSGOps'])
2661    #expand ops to include inversions if any
2662    if SGData['SGInv']:
2663        for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
2664            SGOps.append([-op[0],-op[1]%1.])
2665            SSGOps.append([-sop[0],-sop[1]%1.])
2666    #build set of sym ops around special position       
2667    SSop = []
2668    Sop = []
2669    Sdtau = []
2670    for iop,Op in enumerate(SGOps):         
2671        nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
2672        if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
2673            SSop.append(SSGOps[iop])
2674            Sop.append(SGOps[iop])
2675            ssopinv = nl.inv(SSGOps[iop][0])
2676            mst = ssopinv[3][:3]
2677            epsinv = ssopinv[3][3]
2678            Sdtau.append(np.sum(mst*(XYZ-SGOps[iop][1])-epsinv*SSGOps[iop][1][3]))
2679    SdIndx = np.argsort(np.array(Sdtau))     # just to do in sensible order
2680    if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
2681    #setup displacement arrays
2682    tau = np.linspace(-1,1,49,True)
2683    #make modulation arrays - one parameter at a time
2684    if Stype == 'Sfrac':
2685        CSI,dF,dFTP = DoFrac()
2686    elif Stype == 'Spos':
2687        CSI,dF,dFTP = DoXYZ()
2688        CSI[0] = orderParms(CSI[0])
2689    elif Stype == 'Sadp':
2690        CSI,dF,dFTP = DoUij()
2691        CSI[0] = orderParms(CSI[0]) 
2692    elif Stype == 'Smag':
2693        CSI,dF,dFTP = DoMag()
2694       
2695    if debug:
2696        return CSI,dF,dFTP
2697    else:
2698        return CSI
2699   
2700def MustrainNames(SGData):
2701    'Needs a doc string'
2702    laue = SGData['SGLaue']
2703    uniq = SGData['SGUniq']
2704    if laue in ['m3','m3m']:
2705        return ['S400','S220']
2706    elif laue in ['6/m','6/mmm','3m1']:
2707        return ['S400','S004','S202']
2708    elif laue in ['31m','3']:
2709        return ['S400','S004','S202','S211']
2710    elif laue in ['3R','3mR']:
2711        return ['S400','S220','S310','S211']
2712    elif laue in ['4/m','4/mmm']:
2713        return ['S400','S004','S220','S022']
2714    elif laue in ['mmm']:
2715        return ['S400','S040','S004','S220','S202','S022']
2716    elif laue in ['2/m']:
2717        SHKL = ['S400','S040','S004','S220','S202','S022']
2718        if uniq == 'a':
2719            SHKL += ['S013','S031','S211']
2720        elif uniq == 'b':
2721            SHKL += ['S301','S103','S121']
2722        elif uniq == 'c':
2723            SHKL += ['S130','S310','S112']
2724        return SHKL
2725    else:
2726        SHKL = ['S400','S040','S004','S220','S202','S022']
2727        SHKL += ['S310','S103','S031','S130','S301','S013']
2728        SHKL += ['S211','S121','S112']
2729        return SHKL
2730       
2731def HStrainVals(HSvals,SGData):
2732    laue = SGData['SGLaue']
2733    uniq = SGData['SGUniq']
2734    DIJ = np.zeros(6)
2735    if laue in ['m3','m3m']:
2736        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
2737    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2738        DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
2739    elif laue in ['3R','3mR']:
2740        DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
2741    elif laue in ['4/m','4/mmm']:
2742        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
2743    elif laue in ['mmm']:
2744        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2745    elif laue in ['2/m']:
2746        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2747        if uniq == 'a':
2748            DIJ[5] = HSvals[3]
2749        elif uniq == 'b':
2750            DIJ[4] = HSvals[3]
2751        elif uniq == 'c':
2752            DIJ[3] = HSvals[3]
2753    else:
2754        DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
2755    return DIJ
2756
2757def HStrainNames(SGData):
2758    'Needs a doc string'
2759    laue = SGData['SGLaue']
2760    uniq = SGData['SGUniq']
2761    if laue in ['m3','m3m']:
2762        return ['D11','eA']         #add cubic strain term
2763    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2764        return ['D11','D33']
2765    elif laue in ['3R','3mR']:
2766        return ['D11','D12']
2767    elif laue in ['4/m','4/mmm']:
2768        return ['D11','D33']
2769    elif laue in ['mmm']:
2770        return ['D11','D22','D33']
2771    elif laue in ['2/m']:
2772        Dij = ['D11','D22','D33']
2773        if uniq == 'a':
2774            Dij += ['D23']
2775        elif uniq == 'b':
2776            Dij += ['D13']
2777        elif uniq == 'c':
2778            Dij += ['D12']
2779        return Dij
2780    else:
2781        Dij = ['D11','D22','D33','D12','D13','D23']
2782        return Dij
2783   
2784def MustrainCoeff(HKL,SGData):
2785    'Needs a doc string'
2786    #NB: order of terms is the same as returned by MustrainNames
2787    laue = SGData['SGLaue']
2788    uniq = SGData['SGUniq']
2789    h,k,l = HKL
2790    Strm = []
2791    if laue in ['m3','m3m']:
2792        Strm.append(h**4+k**4+l**4)
2793        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
2794    elif laue in ['6/m','6/mmm','3m1']:
2795        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
2796        Strm.append(l**4)
2797        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
2798    elif laue in ['31m','3']:
2799        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
2800        Strm.append(l**4)
2801        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
2802        Strm.append(4.0*h*k*l*(h+k))
2803    elif laue in ['3R','3mR']:
2804        Strm.append(h**4+k**4+l**4)
2805        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
2806        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
2807        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
2808    elif laue in ['4/m','4/mmm']:
2809        Strm.append(h**4+k**4)
2810        Strm.append(l**4)
2811        Strm.append(3.0*(h*k)**2)
2812        Strm.append(3.0*((h*l)**2+(k*l)**2))
2813    elif laue in ['mmm']:
2814        Strm.append(h**4)
2815        Strm.append(k**4)
2816        Strm.append(l**4)
2817        Strm.append(3.0*(h*k)**2)
2818        Strm.append(3.0*(h*l)**2)
2819        Strm.append(3.0*(k*l)**2)
2820    elif laue in ['2/m']:
2821        Strm.append(h**4)
2822        Strm.append(k**4)
2823        Strm.append(l**4)
2824        Strm.append(3.0*(h*k)**2)
2825        Strm.append(3.0*(h*l)**2)
2826        Strm.append(3.0*(k*l)**2)
2827        if uniq == 'a':
2828            Strm.append(2.0*k*l**3)
2829            Strm.append(2.0*l*k**3)
2830            Strm.append(4.0*k*l*h**2)
2831        elif uniq == 'b':
2832            Strm.append(2.0*l*h**3)
2833            Strm.append(2.0*h*l**3)
2834            Strm.append(4.0*h*l*k**2)
2835        elif uniq == 'c':
2836            Strm.append(2.0*h*k**3)
2837            Strm.append(2.0*k*h**3)
2838            Strm.append(4.0*h*k*l**2)
2839    else:
2840        Strm.append(h**4)
2841        Strm.append(k**4)
2842        Strm.append(l**4)
2843        Strm.append(3.0*(h*k)**2)
2844        Strm.append(3.0*(h*l)**2)
2845        Strm.append(3.0*(k*l)**2)
2846        Strm.append(2.0*k*h**3)
2847        Strm.append(2.0*h*l**3)
2848        Strm.append(2.0*l*k**3)
2849        Strm.append(2.0*h*k**3)
2850        Strm.append(2.0*l*h**3)
2851        Strm.append(2.0*k*l**3)
2852        Strm.append(4.0*k*l*h**2)
2853        Strm.append(4.0*h*l*k**2)
2854        Strm.append(4.0*k*h*l**2)
2855    return Strm
2856
2857def MuShklMean(SGData,Amat,Shkl):
2858   
2859    def genMustrain(xyz,Shkl):
2860        uvw = np.inner(Amat.T,xyz)
2861        Strm = np.array(MustrainCoeff(uvw,SGData))
2862        Sum = np.sum(np.multiply(Shkl,Strm))
2863        Sum = np.where(Sum > 0.01,Sum,0.01)
2864        Sum = np.sqrt(Sum)
2865        return Sum*xyz
2866       
2867    PHI = np.linspace(0.,360.,30,True)
2868    PSI = np.linspace(0.,180.,30,True)
2869    X = np.outer(npcosd(PHI),npsind(PSI))
2870    Y = np.outer(npsind(PHI),npsind(PSI))
2871    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
2872    XYZ = np.dstack((X,Y,Z))
2873    XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
2874    return np.sqrt(np.sum(XYZ**2)/900.)
2875   
2876   
2877def Muiso2Shkl(muiso,SGData,cell):
2878    "this is to convert isotropic mustrain to generalized Shkls"
2879    import GSASIIlattice as G2lat
2880    A = G2lat.cell2AB(cell)[0]
2881   
2882    def minMus(Shkl,muiso,H,SGData,A):
2883        U = np.inner(A.T,H)
2884        S = np.array(MustrainCoeff(U,SGData))
2885        Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,np.newaxis]),axis=0))
2886        rad = np.sqrt(np.sum((Sum[:,np.newaxis]*H)**2,axis=1))
2887        return (muiso-rad)**2
2888       
2889    laue = SGData['SGLaue']
2890    PHI = np.linspace(0.,360.,60,True)
2891    PSI = np.linspace(0.,180.,60,True)
2892    X = np.outer(npsind(PHI),npsind(PSI))
2893    Y = np.outer(npcosd(PHI),npsind(PSI))
2894    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
2895    HKL = np.dstack((X,Y,Z))
2896    if laue in ['m3','m3m']:
2897        S0 = [1000.,1000.]
2898    elif laue in ['6/m','6/mmm','3m1']:
2899        S0 = [1000.,1000.,1000.]
2900    elif laue in ['31m','3']:
2901        S0 = [1000.,1000.,1000.,1000.]
2902    elif laue in ['3R','3mR']:
2903        S0 = [1000.,1000.,1000.,1000.]
2904    elif laue in ['4/m','4/mmm']:
2905        S0 = [1000.,1000.,1000.,1000.]
2906    elif laue in ['mmm']:
2907        S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
2908    elif laue in ['2/m']:
2909        S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
2910    else:
2911        S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000., 
2912            1000.,1000.,0.,0.,0.]
2913    S0 = np.array(S0)
2914    HKL = np.reshape(HKL,(-1,3))
2915    result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
2916    return result[0]
2917       
2918def PackRot(SGOps):
2919    IRT = []
2920    for ops in SGOps:
2921        M = ops[0]
2922        irt = 0
2923        for j in range(2,-1,-1):
2924            for k in range(2,-1,-1):
2925                irt *= 3
2926                irt += M[k][j]
2927        IRT.append(int(irt))
2928    return IRT
2929       
2930def SytSym(XYZ,SGData):
2931    '''
2932    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
2933
2934    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
2935    :param SGData: from SpcGroup
2936    :Returns: a two element tuple:
2937
2938     * The 1st element is a code for the site symmetry (see GetKNsym)
2939     * The 2nd element is the site multiplicity
2940
2941    '''
2942    Mult = 1
2943    Isym = 0
2944    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
2945        Isym = 1073741824
2946    Jdup = 0
2947    Ndup = 0
2948    dupDir = {}
2949    inv = SGData['SGInv']+1
2950    icen = SGData['SGCen']
2951    if SGData['SGFixed']:       #already in list of operators
2952        inv = 1
2953    Xeqv = GenAtom(XYZ,SGData,True)
2954    IRT = PackRot(SGData['SGOps'])
2955    L = -1
2956    for ic,cen in enumerate(icen):
2957        for invers in range(int(inv)):
2958            for io,ops in enumerate(SGData['SGOps']):
2959                irtx = (1-2*invers)*IRT[io]
2960                L += 1
2961                if not Xeqv[L][1]:
2962                    Ndup = io
2963                    Jdup += 1
2964                    jx = GetOprPtrName(str(irtx))   #[KN table no,op name,KNsym ptr]
2965                    if jx[2] < 39:
2966                        px = GetOprName(str(irtx))
2967                        if px != '6643':    #skip Iden
2968                            dupDir[px] = io
2969                        Isym += 2**(jx[2]-1)
2970    if Isym == 1073741824: Isym = 0
2971    Mult = len(SGData['SGOps'])*len(SGData['SGCen'])*(int(SGData['SGInv'])+1)//Jdup
2972         
2973    return GetKNsym(str(Isym)),Mult,Ndup,dupDir
2974   
2975def ElemPosition(SGData):
2976    ''' Under development.
2977    Object here is to return a list of symmetry element types and locations suitable
2978    for say drawing them.
2979    So far I have the element type... getting all possible locations without lookup may be impossible!
2980    '''
2981    Inv = SGData['SGInv']
2982    eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
2983    # get operators & expand if centrosymmetric
2984    Ops = SGData['SGOps']
2985    opM = np.array([op[0].T for op in Ops])
2986    opT = np.array([op[1] for op in Ops])
2987    if Inv:
2988        opM = np.concatenate((opM,-opM))
2989        opT = np.concatenate((opT,-opT))
2990    opMT = list(zip(opM,opT))
2991    for M,T in opMT[1:]:        #skip I
2992        Dt = int(nl.det(M))
2993        Tr = int(np.trace(M))
2994        Dt = -(Dt-1)//2
2995        Es = eleSym[Tr][Dt]
2996        if Dt:              #rotation-inversion
2997            I = np.eye(3)
2998            if Tr == 1:     #mirrors/glides
2999                if np.any(T):       #glide
3000                    M2 = np.inner(M,M)
3001                    MT = np.inner(M,T)+T
3002                    print ('glide',Es,MT)
3003                    print (M2)
3004                else:               #mirror
3005                    print ('mirror',Es,T)
3006                    print (I-M)
3007                X = [-1,-1,-1]
3008            elif Tr == -3:  # pure inversion
3009                X = np.inner(nl.inv(I-M),T)
3010                print ('inversion',Es,X)
3011            else:           #other rotation-inversion
3012                M2 = np.inner(M,M)
3013                MT = np.inner(M,T)+T
3014                print ('rot-inv',Es,MT)
3015                print (M2)
3016                X = [-1,-1,-1]
3017        else:               #rotations
3018            print ('rotation',Es)
3019            X = [-1,-1,-1]
3020        #SymElements.append([Es,X])
3021       
3022    return #SymElements
3023   
3024def ApplyStringOps(A,SGData,X,Uij=[]):
3025    'Needs a doc string'
3026    SGOps = SGData['SGOps']
3027    SGCen = SGData['SGCen']
3028    Ax = A.split('+')
3029    Ax[0] = int(Ax[0])
3030    iC = 0
3031    if Ax[0] < 0:
3032        iC = 1
3033    Ax[0] = abs(Ax[0])
3034    nA = Ax[0]%100-1
3035    cA = Ax[0]//100
3036    Cen = SGCen[cA]
3037    M,T = SGOps[nA]
3038    if len(Ax)>1:
3039        cellA = Ax[1].split(',')
3040        cellA = np.array([int(a) for a in cellA])
3041    else:
3042        cellA = np.zeros(3)
3043    newX = Cen+(1-2*iC)*(np.inner(M,X).T+T)+cellA
3044    if len(Uij):
3045        U = Uij2U(Uij)
3046        U = np.inner(M,np.inner(U,M).T)
3047        newUij = U2Uij(U)
3048        return [newX,newUij]
3049    else:
3050        return newX
3051       
3052def ApplyStringOpsMom(A,SGData,Mom):
3053    'Needs a doc string'
3054    SGOps = SGData['SGOps']
3055    Ax = A.split('+')
3056    Ax[0] = int(Ax[0])
3057    Ax[0] = abs(Ax[0])
3058    nA = Ax[0]%100-1
3059    M,T = SGOps[nA]
3060    if SGData['SGGray']:
3061        newMom = -(np.inner(Mom,M).T)*nl.det(M)
3062    else:
3063        newMom = -(np.inner(Mom,M).T)*SGData['SpnFlp'][nA-1]*nl.det(M)
3064    return newMom
3065       
3066def StringOpsProd(A,B,SGData):
3067    """
3068    Find A*B where A & B are in strings '-' + '100*c+n' + '+ijk'
3069    where '-' indicates inversion, c(>0) is the cell centering operator,
3070    n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3071    Should return resultant string - C. SGData - dictionary using entries:
3072
3073       *  'SGCen': cell centering vectors [0,0,0] at least
3074       *  'SGOps': symmetry operations as [M,T] so that M*x+T = x'
3075
3076    """
3077    SGOps = SGData['SGOps']
3078    SGCen = SGData['SGCen']
3079    #1st split out the cell translation part & work on the operator parts
3080    Ax = A.split('+'); Bx = B.split('+')
3081    Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
3082    iC = 0
3083    if Ax[0]*Bx[0] < 0:
3084        iC = 1
3085    Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
3086    nA = Ax[0]%100-1;  nB = Bx[0]%100-1
3087    cA = Ax[0]//100;  cB = Bx[0]//100
3088    Cen = (SGCen[cA]+SGCen[cB])%1.0
3089    cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
3090    Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
3091    Mc = np.inner(Ma,Mb.T)
3092#    print Ma,Mb,Mc
3093    Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
3094#    print Ta,Tb,Tc
3095#    print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
3096    nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
3097    #now the cell translation part
3098    if len(Ax)>1:
3099        cellA = Ax[1].split(',')
3100        cellA = [int(a) for a in cellA]
3101    else:
3102        cellA = [0,0,0]
3103    if len(Bx)>1:
3104        cellB = Bx[1].split(',')
3105        cellB = [int(b) for b in cellB]
3106    else:
3107        cellB = [0,0,0]
3108    cellC = np.add(cellA,cellB)
3109    C = str(((nC+1)+(100*cC))*(1-2*iC))+'+'+ \
3110        str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
3111    return C
3112           
3113def U2Uij(U):
3114    #returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
3115    return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
3116   
3117def Uij2U(Uij):
3118    #returns the thermal motion tensor U from Uij as numpy array
3119    return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
3120
3121def StandardizeSpcName(spcgroup):
3122    '''Accept a spacegroup name where spaces may have not been used
3123    in the names according to the GSAS convention (spaces between symmetry
3124    for each axis) and return the space group name as used in GSAS
3125    '''
3126    rspc = spcgroup.replace(' ','').upper()
3127    # deal with rhombohedral and hexagonal setting designations
3128    rhomb = ''
3129    if rspc[-1:] == 'R':
3130        rspc = rspc[:-1]
3131        rhomb = ' R'
3132    gray = ''
3133    if "1'" in rspc:
3134        gray = " 1'"
3135        rspc = rspc.replace("1'",'')
3136    elif rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
3137        rspc = rspc[:-1]
3138    else:
3139        rspc = rspc.replace("'",'')
3140    # look for a match in the spacegroup lists
3141    for i in spglist.values():
3142        for spc in i:
3143            if rspc == spc.replace(' ','').upper():
3144                return spc+gray+rhomb
3145    # how about the post-2002 orthorhombic names?
3146    if rspc in sgequiv_2002_orthorhombic:
3147        return sgequiv_2002_orthorhombic[rspc]+gray
3148    else:
3149    # not found
3150        return ''
3151
3152spgbyNum = []
3153'''Space groups indexed by number'''
3154spgbyNum = [None,
3155        'P 1','P -1',                                                   #1-2
3156        'P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
3157        'C 2/m','P 2/c','P 21/c','C 2/c',                               #3-15
3158        'P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21',
3159        'C 2 2 21','C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
3160        'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
3161        'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
3162        'C m m 2','C m c 21','C c c 2',
3163        'A m m 2','A b m 2','A m a 2','A b a 2',
3164        'F m m 2','F d d 2','I m m 2','I b a 2','I m a 2',
3165        'P m m m','P n n n','P c c m','P b a n',
3166        'P m m a','P n n a','P m n a','P c c a','P b a m','P c c n',
3167        'P b c m','P n n m','P m m n','P b c n','P b c a','P n m a',
3168        'C m c m','C m c a','C m m m','C c c m','C m m a','C c c a',
3169        'F m m m', 'F d d d',
3170        'I m m m','I b a m','I b c a','I m m a',                        #16-74
3171        'P 4','P 41','P 42','P 43',
3172        'I 4','I 41',
3173        'P -4','I -4','P 4/m','P 42/m','P 4/n','P 42/n',
3174        'I 4/m','I 41/a',
3175        'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3176        'P 42 21 2','P 43 2 2','P 43 21 2',
3177        'I 4 2 2','I 41 2 2',
3178        'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c','P 4 n c',
3179        'P 42 m c','P 42 b c',
3180        'I 4 m m','I 4 c m','I 41 m d','I 41 c d',
3181        'P -4 2 m','P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2',
3182        'P -4 c 2','P -4 b 2','P -4 n 2',
3183        'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d',
3184        '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',
3185        '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',
3186        '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',
3187        'P 42/n c m',
3188        'I 4/m m m','I 4/m c m','I 41/a m d','I 41/a c d',
3189        'P 3','P 31','P 32','R 3','P -3','R -3',
3190        'P 3 1 2','P 3 2 1','P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1',
3191        'R 3 2',
3192        'P 3 m 1','P 3 1 m','P 3 c 1','P 3 1 c',
3193        'R 3 m','R 3 c',
3194        'P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1',
3195        'R -3 m','R -3 c',                                               #75-167
3196        'P 6','P 61',
3197        'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3198        'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3199        'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3200        '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
3201        'P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m 3','P n 3',
3202        'F m -3','F d -3','I m -3',
3203        'P a -3','I a -3','P 4 3 2','P 42 3 2','F 4 3 2','F 41 3 2',
3204        'I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2','P -4 3 m',
3205        'F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
3206        'P m -3 m','P n -3 n','P m -3 n','P n -3 m',
3207        'F m -3 m','F m -3 c','F d -3 m','F d -3 c',
3208        'I m -3 m','I a -3 d',]                                       #195-230
3209spg2origins = {}
3210''' A dictionary of all spacegroups that have 2nd settings; the value is the
32111st --> 2nd setting transformation vector as X(2nd) = X(1st)-V, nonstandard ones are included.
3212'''
3213spg2origins = {
3214        'P n n n':[-.25,-.25,-.25],
3215        'P b a n':[-.25,-.25,0],'P n c b':[0,-.25,-.25],'P c n a':[-.25,0,-.25],
3216        'P m m n':[-.25,-.25,0],'P n m m':[0,-.25,-.25],'P m n m':[-.25,0,-.25],
3217        'C c c a':[0,-.25,-.25],'C c c b':[-.25,0,-.25],'A b a a':[-.25,0,-.25],
3218        'A c a a':[-.25,-.25,0],'B b c b':[-.25,-.25,0],'B b a b':[0,-.25,-.25],
3219        'F d d d':[-.125,-.125,-.125],
3220        'P 4/n':[-.25,-.25,0],'P 42/n':[-.25,-.25,-.25],'I 41/a':[0,-.25,-.125],
3221        '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],
3222        '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],
3223        'I 41/a m d':[0,.25,-.125],'I 41/a c d':[0,.25,-.125],
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 -3 c':[-.375,-.375,-.375],
3226        'p n 3':[-.25,-.25,-.25],'F d 3':[-.125,-.125,-.125],'P n 3 n':[-.25,-.25,-.25],
3227        'P n 3 m':[-.25,-.25,-.25],'F d 3 m':[-.125,-.125,-.125],'F d - c':[-.375,-.375,-.375]}
3228spglist = {}
3229'''A dictionary of space groups as ordered and named in the pre-2002 International
3230Tables Volume A, except that spaces are used following the GSAS convention to
3231separate the different crystallographic directions.
3232Note that the symmetry codes here will recognize many non-standard space group
3233symbols with different settings. They are ordered by Laue group
3234'''
3235spglist = {
3236    'P1' : ('P 1','P -1',), # 1-2
3237    'C1' : ('C 1','C -1',),
3238    'P2/m': ('P 2','P 21','P m','P a','P c','P n',
3239        '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
3240    'C2/m':('C 2','C m','C c','C n',
3241        'C 2/m','C 2/c','C 2/n',),
3242    'A2/m':('A 2','A m','A a','A n',
3243        'A 2/m','A 2/a','A 2/n',),
3244   'Pmmm':('P 2 2 2',
3245        'P 2 2 21','P 21 2 2','P 2 21 2',
3246        'P 21 21 2','P 2 21 21','P 21 2 21',
3247        'P 21 21 21',
3248        'P m m 2','P 2 m m','P m 2 m',
3249        'P m c 21','P 21 m a','P b 21 m','P m 21 b','P c m 21','P 21 a m',
3250        'P c c 2','P 2 a a','P b 2 b',
3251        'P m a 2','P 2 m b','P c 2 m','P m 2 a','P b m 2','P 2 c m',
3252        'P c a 21','P 21 a b','P c 21 b','P b 21 a','P b c 21','P 21 c a',
3253        'P n c 2','P 2 n a','P b 2 n','P n 2 b','P c n 2','P 2 a n',
3254        'P m n 21','P 21 m n','P n 21 m','P m 21 n','P n m 21','P 21 n m',
3255        'P b a 2','P 2 c b','P c 2 a',
3256        'P n a 21','P 21 n b','P c 21 n','P n 21 a','P b n 21','P 21 c n',
3257        'P n n 2','P 2 n n','P n 2 n',
3258        'P m m m','P n n n',
3259        'P c c m','P m a a','P b m b',
3260        'P b a n','P n c b','P c n a',
3261        'P m m a','P b m m','P m c m','P m a m','P m m b','P c m m',
3262        'P n n a','P b n n','P n c n','P n a n','P n n b','P c n n',
3263        'P m n a','P b m n','P n c m','P m a n','P n m b','P c n m',
3264        'P c c a','P b a a','P b c b','P b a b','P c c b','P c a a',
3265        'P b a m','P m c b','P c m a',
3266        'P c c n','P n a a','P b n b',
3267        'P b c m','P m c a','P b m a','P c m b','P c a m','P m a b',
3268        'P n n m','P m n n','P n m n',
3269        'P m m n','P n m m','P m n m',
3270        'P b c n','P n c a','P b n a','P c n b','P c a n','P n a b',
3271        'P b c a','P c a b',
3272        'P n m a','P b n m','P m c n','P n a m','P m n b','P c m n',
3273        ),
3274    'Cmmm':('C 2 2 21','C 2 2 2','C m m 2',
3275        'C m c 21','C c m 21','C c c 2','C m 2 m','C 2 m m',
3276        'C m 2 a','C 2 m b','C c 2 m','C 2 c m','C c 2 a','C 2 c b',
3277        'C m c m','C c m m','C m c a','C c m b',
3278        'C m m m','C c c m','C m m a','C m m b','C c c a','C c c b',),
3279    'Ammm':('A 21 2 2','A 2 2 2','A 2 m m',
3280        'A 21 m a','A 21 a m','A 2 a a','A m 2 m','A m m 2',
3281        'A b m 2','A c 2 m','A m a 2','A m 2 a','A b a 2','A c 2 a',
3282        'A m m a','A m a m','A b m a','A c a m',
3283        'A m m m','A m a a','A b m m','A c m m','A c a a','A b a a',),
3284    'Bmmm':('B 2 21 2','B 2 2 2','B m 2 m',
3285        'B m 21 b','B b 21 m','B b 2 b','B m m 2','B 2 m m',
3286        'B 2 c m','B m a 2','B 2 m b','B b m 2','B 2 c b','B b a 2',
3287        'B b m m','B m m b','B b c m','B m a b',
3288        'B m m m','B b m b','B m a m','B m c m','B b a b','B b c b',),
3289    'Immm':('I 2 2 2','I 21 21 21',
3290        'I m m 2','I m 2 m','I 2 m m',
3291        'I b a 2','I 2 c b','I c 2 a',
3292        'I m a 2','I 2 m b','I c 2 m','I m 2 a','I b m 2','I 2 c m',
3293        'I m m m','I b a m','I m c b','I c m a',
3294        'I b c a','I c a b',
3295        'I m m a','I b m m ','I m c m','I m a m','I m m b','I c m m',),
3296    'Fmmm':('F 2 2 2','F m m m', 'F d d d',
3297        'F m m 2','F m 2 m','F 2 m m',
3298        'F d d 2','F d 2 d','F 2 d d',),
3299    'P4/mmm':('P 4','P 41','P 42','P 43','P -4','P 4/m','P 42/m','P 4/n','P 42/n',
3300        'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3301        'P 42 21 2','P 43 2 2','P 43 21 2','P 4 m m','P 4 b m','P 42 c m',
3302        'P 42 n m','P 4 c c','P 4 n c','P 42 m c','P 42 b c','P -4 2 m',
3303        'P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2',
3304        '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',
3305        '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',
3306        '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',
3307        'P 42/n c m',),
3308    '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',
3309        'I 4 c m','I 41 m d','I 41 c d',
3310        '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',
3311        'I 41/a m d','I 41/a c d'),
3312    'R3-H':('R 3','R -3','R 3 2','R 3 m','R 3 c','R -3 m','R -3 c',),
3313    'P6/mmm': ('P 3','P 31','P 32','P -3','P 3 1 2','P 3 2 1','P 31 1 2',
3314        'P 31 2 1','P 32 1 2','P 32 2 1', 'P 3 m 1','P 3 1 m','P 3 c 1',
3315        'P 3 1 c','P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1','P 6','P 61',
3316        'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3317        'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3318        'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3319        'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',),
3320    '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',
3321        'P 4 3 2','P 42 3 2','P 43 3 2','P 41 3 2','P -4 3 m','P -4 3 n',
3322        'P m 3 m','P m -3 m','P n 3 n','P n -3 n',
3323        'P m 3 n','P m -3 n','P n 3 m','P n -3 m',),
3324    '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',
3325        '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'),
3326    'Fm3m':('F 2 3','F m 3','F m -3','F d 3','F d -3',
3327        'F 4 3 2','F 41 3 2','F -4 3 m','F -4 3 c',
3328        'F m 3 m','F m -3 m','F m 3 c','F m -3 c',
3329        'F d 3 m','F d -3 m','F d 3 c','F d -3 c',),
3330}
3331sgequiv_2002_orthorhombic = {}
3332''' A dictionary of orthorhombic space groups that were renamed in the 2002 Volume A,
3333 along with the pre-2002 name. The e designates a double glide-plane
3334'''
3335sgequiv_2002_orthorhombic = {
3336        'AEM2':'A b m 2','B2EM':'B 2 c m','CM2E':'C m 2 a',
3337        'AE2M':'A c 2 m','BME2':'B m a 2','C2ME':'C 2 m b',
3338        'AEA2':'A b a 2','B2EB':'B 2 c b','CC2E':'C c 2 a',
3339        'AE2A':'A c 2 a','BBE2':'B b a 2','C2CE':'C 2 c b',
3340        'CMCE':'C m c a','AEMA':'A b m a','BBEM':'B b c m',
3341        'BMEB':'B m a b','CCME':'C c m b','AEAM':'A c a m',
3342        'CMME':'C m m a','AEMM':'A b m m','BMEM':'B m c m',
3343        'CCCE':'C c c a','AEAA':'A b a a','BBEB':'B b c b'}
3344
3345#'A few non-standard space groups for test use'
3346nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h', 
3347                      'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
3348                      'R 3 c r','R -3 c r','R -3 m r',),
3349
3350#Use the space groups types in this order to list the symbols in the
3351#order they are listed in the International Tables, vol. A'''
3352symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 
3353               'trigonal', 'hexagonal', 'cubic')
3354
3355# self-test materials follow. Requires files in directory testinp
3356selftestlist = []
3357'''Defines a list of self-tests'''
3358selftestquiet = True
3359def _ReportTest():
3360    'Report name and doc string of current routine when ``selftestquiet`` is False'
3361    if not selftestquiet:
3362        import inspect
3363        caller = inspect.stack()[1][3]
3364        doc = eval(caller).__doc__
3365        if doc is not None:
3366            print('testing '+__file__+' with '+caller+' ('+doc+')')
3367        else:
3368            print('testing '+__file__()+" with "+caller)
3369def test0():
3370    '''self-test #0: exercise MoveToUnitCell'''
3371    _ReportTest()
3372    msg = "MoveToUnitCell failed"
3373    assert (MoveToUnitCell([1,2,3])[0] == [0,0,0]).all, msg
3374    assert (MoveToUnitCell([2,-1,-2])[0] == [0,0,0]).all, msg
3375    assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9)[0] < 1e-6, msg
3376    assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1)[0] < 1e-6, msg
3377selftestlist.append(test0)
3378
3379def test1():
3380    '''self-test #1: SpcGroup against previous results'''
3381    #'''self-test #1: SpcGroup and SGPrint against previous results'''
3382    _ReportTest()
3383    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3384    if ospath.exists(testdir):
3385        if testdir not in sys.path: sys.path.insert(0,testdir)
3386    import spctestinp
3387    def CompareSpcGroup(spc, referr, refdict, reflist): 
3388        'Compare output from GSASIIspc.SpcGroup with results from a previous run'
3389        # if an error is reported, the dictionary can be ignored
3390        msg0 = "CompareSpcGroup failed on space group %s" % spc
3391        result = SpcGroup(spc)
3392        if result[0] == referr and referr > 0: return True
3393#        #print result[1]['SpGrp']
3394        #msg = msg0 + " in list lengths"
3395        #assert len(keys) == len(refdict.keys()), msg
3396        for key in refdict.keys():
3397            if key == 'SGOps' or  key == 'SGCen':
3398                msg = msg0 + (" in key %s length" % key)
3399                assert len(refdict[key]) == len(result[1][key]), msg
3400                for i in range(len(refdict[key])):
3401                    msg = msg0 + (" in key %s level 0" % key)
3402                    assert np.allclose(result[1][key][i][0],refdict[key][i][0]), msg
3403                    msg = msg0 + (" in key %s level 1" % key)
3404                    assert np.allclose(result[1][key][i][1],refdict[key][i][1]), msg
3405            else:
3406                msg = msg0 + (" in key %s" % key)
3407                assert result[1][key] == refdict[key], msg
3408        msg = msg0 + (" in key %s reflist" % key)
3409        #for (l1,l2) in zip(reflist, SGPrint(result[1])):
3410        #    assert l2.replace('\t','').replace(' ','') == l1.replace(' ',''), 'SGPrint ' +msg
3411        # for now disable SGPrint testing, output has changed
3412        #assert reflist == SGPrint(result[1]), 'SGPrint ' +msg
3413    for spc in spctestinp.SGdat:
3414        CompareSpcGroup(spc, 0, spctestinp.SGdat[spc], spctestinp.SGlist[spc] )
3415selftestlist.append(test1)
3416
3417def test2():
3418    '''self-test #2: SpcGroup against cctbx (sgtbx) computations'''
3419    _ReportTest()
3420    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3421    if ospath.exists(testdir):
3422        if testdir not in sys.path: sys.path.insert(0,testdir)
3423    import sgtbxtestinp
3424    def CompareWcctbx(spcname, cctbx_in, debug=0):
3425        'Compare output from GSASIIspc.SpcGroup with results from cctbx.sgtbx'
3426        cctbx = cctbx_in[:] # make copy so we don't delete from the original
3427        spc = (SpcGroup(spcname))[1]
3428        if debug: print (spc['SpGrp'])
3429        if debug: print (spc['SGCen'])
3430        latticetype = spcname.strip().upper()[0]
3431        # lattice type of R implies Hexagonal centering", fix the rhombohedral settings
3432        if latticetype == "R" and len(spc['SGCen']) == 1: latticetype = 'P'
3433        assert latticetype == spc['SGLatt'], "Failed: %s does not match Lattice: %s" % (spcname, spc['SGLatt'])
3434        onebar = [1]
3435        if spc['SGInv']: onebar.append(-1)
3436        for (op,off) in spc['SGOps']:
3437            for inv in onebar:
3438                for cen in spc['SGCen']:
3439                    noff = off + cen
3440                    noff = MoveToUnitCell(noff)[0]
3441                    mult = tuple((op*inv).ravel().tolist())
3442                    if debug: print ("\n%s: %s + %s" % (spcname,mult,noff))
3443                    for refop in cctbx:
3444                        if debug: print (refop)
3445                        # check the transform
3446                        if refop[:9] != mult: continue
3447                        if debug: print ("mult match")
3448                        # check the translation
3449                        reftrans = list(refop[-3:])
3450                        reftrans = MoveToUnitCell(reftrans)[0]
3451                        if all(abs(noff - reftrans) < 1.e-5):
3452                            cctbx.remove(refop)
3453                            break
3454                    else:
3455                        assert False, "failed on %s:\n\t %s + %s" % (spcname,mult,noff)
3456    for key in sgtbxtestinp.sgtbx:
3457        CompareWcctbx(key, sgtbxtestinp.sgtbx[key])
3458selftestlist.append(test2)
3459
3460def test3(): 
3461    '''self-test #3: exercise SytSym (includes GetOprPtrName, GenAtom, GetKNsym)
3462     for selected space groups against info in IT Volume A '''
3463    _ReportTest()
3464    def ExerciseSiteSym (spc, crdlist):
3465        'compare site symmetries and multiplicities for a specified space group'
3466        msg = "failed on site sym test for %s" % spc
3467        (E,S) = SpcGroup(spc)
3468        assert not E, msg
3469        for t in crdlist:
3470            symb, m, n, od = SytSym(t[0],S)
3471            if symb.strip() != t[2].strip() or m != t[1]:
3472                print (spc,t[0],m,n,symb,t[2],od)
3473            assert m == t[1]
3474            #assert symb.strip() == t[2].strip()
3475
3476    ExerciseSiteSym('p 1',[
3477            ((0.13,0.22,0.31),1,'1'),
3478            ((0,0,0),1,'1'),
3479            ])
3480    ExerciseSiteSym('p -1',[
3481            ((0.13,0.22,0.31),2,'1'),
3482            ((0,0.5,0),1,'-1'),
3483            ])
3484    ExerciseSiteSym('C 2/c',[
3485            ((0.13,0.22,0.31),8,'1'),
3486            ((0.0,.31,0.25),4,'2(y)'),
3487            ((0.25,.25,0.5),4,'-1'),
3488            ((0,0.5,0),4,'-1'),
3489            ])
3490    ExerciseSiteSym('p 2 2 2',[
3491            ((0.13,0.22,0.31),4,'1'),
3492            ((0,0.5,.31),2,'2(z)'),
3493            ((0.5,.31,0.5),2,'2(y)'),
3494            ((.11,0,0),2,'2(x)'),
3495            ((0,0.5,0),1,'222'),
3496            ])
3497    ExerciseSiteSym('p 4/n',[
3498            ((0.13,0.22,0.31),8,'1'),
3499            ((0.25,0.75,.31),4,'2(z)'),
3500            ((0.5,0.5,0.5),4,'-1'),
3501            ((0,0.5,0),4,'-1'),