source: trunk/GSASIIspc.py @ 3197

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

Rearrange SingleStringDialog? window - put '?' at end of text box
modify space group stuff & cif importer to read any magnetic cif file.

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