source: trunk/GSASIIspc.py @ 3215

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

mag incommensurate stuff

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