source: trunk/GSASIIspc.py @ 3220

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

minor fixes to commensurate magnetic cif import & display

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