source: trunk/GSASIIspc.py @ 3219

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

fix import mag commensurate structures - all drawn correctly now
modulated mags not correct yet

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