source: trunk/GSASIIspc.py @ 3211

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

fix integer divide in wave display
some fixes to mcif import

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