source: trunk/GSASIIspc.py @ 3212

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

fix problem of missing rd.MPhase in phase imports
work on display of gray phases & other magnetic symmetry info.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 147.5 KB
Line 
1# -*- coding: utf-8 -*-
2"""
3*GSASIIspc: Space group module*
4-------------------------------
5
6Space group interpretation routines. Note that space group information is
7stored in a :ref:`Space Group (SGData)<SGData_table>` object.
8
9"""
10########### SVN repository information ###################
11# $Date: 2018-01-04 16:29:00 +0000 (Thu, 04 Jan 2018) $
12# $Author: vondreele $
13# $Revision: 3212 $
14# $URL: trunk/GSASIIspc.py $
15# $Id: GSASIIspc.py 3212 2018-01-04 16:29:00Z 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: 3212 $")
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    SSgSpc = SSGData['SSpGrp']
1575    if SGData.get('SGGray',False):
1576        SSgSpc = SSgSpc.replace('('," 1'(")
1577    SSGText.append(' Superspace Group: '+SSgSpc)
1578    CentStr = 'centrosymmetric'
1579    if not SGData['SGInv']:
1580        CentStr = 'non'+CentStr
1581    if SGData['SGLatt'] in 'ABCIFR':
1582        SSGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
1583    else:
1584        SSGText.append(' The superlattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())
1585    SSGText.append(' The Laue symmetry is '+SGData['SGLaue'])
1586    SGptGp = SGData['SGPtGrp']
1587    if SGData.get('SGGray',False):
1588        SGptGp += "1'"
1589    SSGText.append(' The superlattice point group is '+SGptGp+', '+''.join([str(i) for i in SSGData['SSGKl']]))
1590    SSGText.append(' The number of superspace group generators is '+str(len(SGData['SSGKl'])))
1591    SSGText.append(' Multiplicity of a general site is '+str(Mult))
1592    if SGData['SGUniq'] in ['a','b','c']:
1593        SSGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
1594    if SGData['SGInv']:
1595        SSGText.append(' The inversion center is located at 0,0,0')
1596    if SGData['SGPolax']:
1597        SSGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
1598    SSGText.append(' ')
1599    if len(SSGData['SSGCen']) > 1:
1600        SSGText.append(' The equivalent positions are:')
1601        SSGText.append(' ('+SSLatt2text(SSGData['SSGCen'])+')+\n')
1602    else:
1603        SSGText.append(' The equivalent positions are:\n')
1604    SSGTable = []
1605    for i,Opr in enumerate(SSGData['SSGOps']):
1606        SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr)))
1607    if SGData.get('SGGray',False):
1608        SSGTable.append("     for 1'")
1609        for i,Opr in enumerate(SSGData['SSGOps']):
1610            Opr2 = [Opr[0],Opr[1]+np.array([0,0,0,.5])]
1611            SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr2)))
1612    return SSGText,SSGTable
1613   
1614def SSGModCheck(Vec,modSymb,newMod=True):
1615    ''' Checks modulation vector compatibility with supersymmetry space group symbol.
1616    if newMod: Superspace group symbol takes precidence & the vector will be modified accordingly
1617    '''
1618    Fracs = {'1/2':0.5,'1/3':1./3,'1':1.0,'0':0.,'a':0.,'b':0.,'g':0.}
1619    modQ = [Fracs[mod] for mod in modSymb]
1620    if newMod:
1621        newVec = [0.1 if (vec == 0.0 and mod in ['a','b','g']) else vec for [vec,mod] in zip(Vec,modSymb)]
1622        return [Q if mod not in ['a','b','g'] and vec != Q else vec for [vec,mod,Q] in zip(newVec,modSymb,modQ)],  \
1623            [True if mod in ['a','b','g'] else False for mod in modSymb]
1624    else:
1625        return Vec,[True if mod in ['a','b','g'] else False for mod in modSymb]
1626
1627def SSMT2text(Opr):
1628    "From superspace group matrix/translation operator returns text version"
1629    XYZS = ('x','y','z','t')    #Stokes, Campbell & van Smaalen notation
1630    TRA = ('   ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
1631    Fld = ''
1632    M,T = Opr
1633    for j in range(4):
1634        IJ = ''
1635        for k in range(4):
1636            txt = str(int(round(M[j][k])))
1637            txt = txt.replace('1',XYZS[k]).replace('0','')
1638            if '2' in txt:
1639                txt += XYZS[k]
1640            if IJ and M[j][k] > 0:
1641                IJ += '+'+txt
1642            else:
1643                IJ += txt
1644        IK = int(round(T[j]*12))%12
1645        if IK:
1646            if not IJ:
1647                break
1648            if IJ[0] == '-':
1649                Fld += (TRA[IK]+IJ).rjust(8)
1650            else:
1651                Fld += (TRA[IK]+'+'+IJ).rjust(8)
1652        else:
1653            Fld += IJ.rjust(8)
1654        if j != 3: Fld += ', '
1655    return Fld
1656   
1657def SSLatt2text(SSGCen):
1658    "Lattice centering vectors to text"
1659    lattTxt = ''
1660    lattDir = {4:'1/3',6:'1/2',8:'2/3',0:'0'}
1661    for vec in SSGCen:
1662        lattTxt += ' '
1663        for item in vec:
1664            lattTxt += '%s,'%(lattDir[int(item*12)])
1665        lattTxt = lattTxt.rstrip(',')
1666        lattTxt += ';'
1667    lattTxt = lattTxt.rstrip(';').lstrip(' ')
1668    return lattTxt
1669       
1670def SSpaceGroup(SGSymbol,SSymbol):
1671    '''
1672    Print the output of SSpcGroup in a nicely formatted way.
1673
1674    :param SGSymbol: space group symbol with spaces between axial fields.
1675    :param SSymbol: superspace group symbol extension (string).
1676    :returns: nothing
1677    '''
1678
1679    E,A = SpcGroup(SGSymbol)
1680    if E > 0:
1681        print (SGErrors(E))
1682        return
1683    E,B = SSpcGroup(A,SSymbol)   
1684    if E > 0:
1685        print (E)
1686        return
1687    for l in SSGPrint(B):
1688        print (l)
1689       
1690def SGProd(OpA,OpB):
1691    '''
1692    Form space group operator product. OpA & OpB are [M,V] pairs;
1693        both must be of same dimension (3 or 4). Returns [M,V] pair
1694    '''
1695    A,U = OpA
1696    B,V = OpB
1697    M = np.inner(B,A.T)
1698    W = np.inner(B,U)+V
1699    return M,W
1700       
1701def MoveToUnitCell(xyz):
1702    '''
1703    Translates a set of coordinates so that all values are >=0 and < 1
1704
1705    :param xyz: a list or numpy array of fractional coordinates
1706    :returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
1707    '''
1708    XYZ = (np.array(xyz)+10.)%1.
1709    cell = np.asarray(np.rint(xyz-XYZ),dtype=np.int32)
1710    return XYZ,cell
1711       
1712def Opposite(XYZ,toler=0.0002):
1713    '''
1714    Gives opposite corner, edge or face of unit cell for position within tolerance.
1715        Result may be just outside the cell within tolerance
1716
1717    :param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
1718    :param toler: unit cell fraction tolerance making opposite
1719    :returns:
1720        XYZ: dict of opposite positions; key=unit cell & always contains XYZ
1721    '''
1722    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]]
1723    TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
1724    perm = TB*perm3
1725    cperm = ['%d,%d,%d'%(i,j,k) for i,j,k in perm]
1726    D = dict(zip(cperm,perm))
1727    new = {}
1728    for key in D:
1729        new[key] = np.array(D[key])+np.array(XYZ)
1730    return new
1731       
1732def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
1733    '''
1734    Generates the equivalent positions for a specified coordinate and space group
1735
1736    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
1737    :param SGData: from :func:`SpcGroup`
1738    :param All: True return all equivalent positions including duplicates;
1739      False return only unique positions
1740    :param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
1741    :param Move: True move generated atom positions to be inside cell
1742      False do not move atoms       
1743    :return: [[XYZEquiv],Idup,[UijEquiv]]
1744
1745      *  [XYZEquiv] is list of equivalent positions (XYZ is first entry)
1746      *  Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
1747      * is centering operator number (1-4) and - is for inversion
1748        Cell = unit cell translations needed to put new positions inside cell
1749        [UijEquiv] - equivalent Uij; absent if no Uij given
1750       
1751    '''
1752    XYZEquiv = []
1753    UijEquiv = []
1754    Idup = []
1755    Cell = []
1756    X = np.array(XYZ)
1757    if Move:
1758        X = MoveToUnitCell(X)[0]
1759    for ic,cen in enumerate(SGData['SGCen']):
1760        C = np.array(cen)
1761        for invers in range(int(SGData['SGInv']+1)):
1762            for io,[M,T] in enumerate(SGData['SGOps']):
1763                idup = ((io+1)+100*ic)*(1-2*invers)
1764                XT = np.inner(M,X)+T
1765                if len(Uij):
1766                    U = Uij2U(Uij)
1767                    U = np.inner(M,np.inner(U,M).T)
1768                    newUij = U2Uij(U)
1769                if invers:
1770                    XT = -XT
1771                XT += C
1772                cell = np.zeros(3,dtype=np.int32)
1773                cellj = np.zeros(3,dtype=np.int32)
1774                if Move:
1775                    newX,cellj = MoveToUnitCell(XT)
1776                else:
1777                    newX = XT
1778                cell += cellj
1779                if All:
1780                    if np.allclose(newX,X,atol=0.0002):
1781                        idup = False
1782                else:
1783                    if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
1784                        idup = False
1785                if All or idup:
1786                    XYZEquiv.append(newX)
1787                    Idup.append(idup)
1788                    Cell.append(cell)
1789                    if len(Uij):
1790                        UijEquiv.append(newUij)                   
1791    if len(Uij):
1792        return list(zip(XYZEquiv,UijEquiv,Idup,Cell))
1793    else:
1794        return list(zip(XYZEquiv,Idup,Cell))
1795       
1796def GenHKL(HKL,SGData):
1797    ''' Generates all equivlent reflections including Friedel pairs
1798    :param HKL:  [h,k,l] must be integral values
1799    :param SGData: space group data obtained from SpcGroup
1800    :returns: array Uniq: equivalent reflections
1801    '''
1802   
1803    Ops = SGData['SGOps']
1804    OpM = np.array([op[0] for op in Ops])
1805    Uniq = np.inner(OpM,HKL)
1806    Uniq = list(Uniq)+list(-1*Uniq)
1807    return np.array(Uniq)
1808
1809def GenHKLf(HKL,SGData):
1810    '''
1811    Uses old GSAS Fortran routine genhkl.for
1812
1813    :param HKL:  [h,k,l] must be integral values for genhkl.for to work
1814    :param SGData: space group data obtained from SpcGroup
1815    :returns: iabsnt,mulp,Uniq,phi
1816
1817     *   iabsnt = True if reflection is forbidden by symmetry
1818     *   mulp = reflection multiplicity including Friedel pairs
1819     *   Uniq = numpy array of equivalent hkl in descending order of h,k,l
1820     *   phi = phase offset for each equivalent h,k,l
1821
1822    '''
1823    hklf = list(HKL)+[0,]       #could be numpy array!
1824    Ops = SGData['SGOps']
1825    OpM = np.array([op[0] for op in Ops],order='F')
1826    OpT = np.array([op[1] for op in Ops])
1827    Cen = np.array([cen for cen in SGData['SGCen']],order='F')
1828   
1829    import pyspg
1830    Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
1831    h,k,l,f = Uniq
1832    Uniq=np.array(list(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq])))
1833    phi = f[:Nuniq]
1834    return iabsnt,mulp,Uniq,phi
1835   
1836def checkSSLaue(HKL,SGData,SSGData):
1837    #Laue check here - Toss HKL if outside unique Laue part
1838    h,k,l,m = HKL
1839    if SGData['SGLaue'] == '2/m':
1840        if SGData['SGUniq'] == 'a':
1841            if 'a' in SSGData['modSymb'] and h == 0 and m < 0:
1842                return False
1843            elif 'b' in SSGData['modSymb'] and k == 0 and l ==0 and m < 0:
1844                return False
1845            else:
1846                return True
1847        elif SGData['SGUniq'] == 'b':
1848            if 'b' in SSGData['modSymb'] and k == 0 and m < 0:
1849                return False
1850            elif 'a' in SSGData['modSymb'] and h == 0 and l ==0 and m < 0:
1851                return False
1852            else:
1853                return True
1854        elif SGData['SGUniq'] == 'c':
1855            if 'g' in SSGData['modSymb'] and l == 0 and m < 0:
1856                return False
1857            elif 'a' in SSGData['modSymb'] and h == 0 and k ==0 and m < 0:
1858                return False
1859            else:
1860                return True
1861    elif SGData['SGLaue'] == 'mmm':
1862        if 'a' in SSGData['modSymb']:
1863            if h == 0 and m < 0:
1864                return False
1865            else:
1866                return True
1867        elif 'b' in SSGData['modSymb']:
1868            if k == 0 and m < 0:
1869                return False
1870            else:
1871                return True
1872        elif 'g' in SSGData['modSymb']:
1873            if l == 0 and m < 0:
1874                return False
1875            else:
1876                return True
1877    else:   #tetragonal, trigonal, hexagonal (& triclinic?)
1878        if l == 0 and m < 0:
1879            return False
1880        else:
1881            return True
1882       
1883   
1884def checkSSextc(HKL,SSGData):
1885    Ops = SSGData['SSGOps']
1886    OpM = np.array([op[0] for op in Ops])
1887    OpT = np.array([op[1] for op in Ops])
1888    HKLS = np.array([HKL,-HKL])     #Freidel's Law
1889    DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,4))
1890    PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
1891    for dhkl,phkl in zip(DHKL,PHKL)[1:]:    #skip identity
1892        if dhkl.any():
1893            continue
1894        else:
1895            if phkl%1.:
1896                return False
1897    return True
1898   
1899################################################################################
1900#### Site symmetry tables
1901################################################################################
1902   
1903OprPtrName = {
1904    '-6643':[   2,' 1bar ', 1],'6479' :[  10,'  2z  ', 2],'-6479':[   9,'  mz  ', 3],
1905    '6481' :[   7,'  my  ', 4],'-6481':[   6,'  2y  ', 5],'6641' :[   4,'  mx  ', 6],
1906    '-6641':[   3,'  2x  ', 7],'6591' :[  28,' m+-0 ', 8],'-6591':[  27,' 2+-0 ', 9],
1907    '6531' :[  25,' m110 ',10],'-6531':[  24,' 2110 ',11],'6537' :[  61,'  4z  ',12],
1908    '-6537':[  62,' -4z  ',13],'975'  :[  68,' 3+++1',14],'6456' :[ 114,'  3z1 ',15],
1909    '-489' :[  73,' 3+-- ',16],'483'  :[  78,' 3-+- ',17],'-969' :[  83,' 3--+ ',18],
1910    '819'  :[  22,' m+0- ',19],'-819' :[  21,' 2+0- ',20],'2431' :[  16,' m0+- ',21],
1911    '-2431':[  15,' 20+- ',22],'-657' :[  19,' m101 ',23],'657'  :[  18,' 2101 ',24],
1912    '1943' :[  48,' -4x  ',25],'-1943':[  47,'  4x  ',26],'-2429':[  13,' m011 ',27],
1913    '2429' :[  12,' 2011 ',28],'639'  :[  55,' -4y  ',29],'-639' :[  54,'  4y  ',30],
1914    '-6484':[ 146,' 2010 ', 4],'6484' :[ 139,' m010 ', 5],'-6668':[ 145,' 2100 ', 6],
1915    '6668' :[ 138,' m100 ', 7],'-6454':[ 148,' 2120 ',18],'6454' :[ 141,' m120 ',19],
1916    '-6638':[ 149,' 2210 ',20],'6638' :[ 142,' m210 ',21],              #search ends here
1917    '2223' :[  68,' 3+++2',39],
1918    '6538' :[ 106,'  6z1 ',40],'-2169':[  83,' 3--+2',41],'2151' :[  73,' 3+--2',42],
1919    '2205' :[  79,'-3-+-2',43],'-2205':[  78,' 3-+-2',44],'489'  :[  74,'-3+--1',45],
1920    '801'  :[  53,'  4y1 ',46],'1945' :[  47,'  4x3 ',47],'-6585':[  62,' -4z3 ',48],
1921    '6585' :[  61,'  4z3 ',49],'6584' :[ 114,'  3z2 ',50],'6666' :[ 106,'  6z5 ',51],
1922    '6643' :[   1,' Iden ',52],'-801' :[  55,' -4y1 ',53],'-1945':[  48,' -4x3 ',54],
1923    '-6666':[ 105,' -6z5 ',55],'-6538':[ 105,' -6z1 ',56],'-2223':[  69,'-3+++2',57],
1924    '-975' :[  69,'-3+++1',58],'-6456':[ 113,' -3z1 ',59],'-483' :[  79,'-3-+-1',60],
1925    '969'  :[  84,'-3--+1',61],'-6584':[ 113,' -3z2 ',62],'2169' :[  84,'-3--+2',63],
1926    '-2151':[  74,'-3+--2',64],'0':[0,' ????',0]
1927    }
1928   
1929OprName = {
1930    '-6643':'   -1   ','6479' :'    2(z)','-6479':'    m(z)',
1931    '6481' :'    m(y)','-6481':'    2(y)','6641' :'    m(x)',
1932    '-6641':'    2(x)','6591' :'  m(+-0)','-6591':'  2(+-0)',
1933    '6531' :' m(110) ','-6531':' 2(110) ','6537' :'  4(001)',
1934    '-6537':' -4(001)','975'  :'  3(111)','6456' :'    3   ',
1935    '-489' :'  3(+--)','483'  :'  3(-+-)','-969' :'  3(--+)',
1936    '819'  :'  m(+0-)','-819' :'  2(+0-)','2431' :'  m(0+-)',
1937    '-2431':'  2(0+-)','-657' :'   m(xz)','657'  :'   2(xz)',
1938    '1943' :' -4(100)','-1943':'  4(100)','-2429':'   m(yz)',
1939    '2429' :'   2(yz)','639'  :' -4(010)','-639' :'  4(010)',
1940    '-6484':' 2(010) ','6484' :' m(010) ','-6668':' 2(100) ',
1941    '6668' :' m(100) ','-6454':' 2(120) ','6454' :' m(120) ',
1942    '-6638':' 2(210) ','6638' :' m(210) '}              #search ends here
1943   
1944                                 
1945KNsym = {
1946    '0'         :'    1   ','1'         :'   -1   ','64'        :'    2(x)','32'        :'    m(x)',
1947    '97'        :'  2/m(x)','16'        :'    2(y)','8'         :'    m(y)','25'        :'  2/m(y)',
1948    '2'         :'    2(z)','4'         :'    m(z)','7'         :'  2/m(z)','134217728' :'   2(yz)',
1949    '67108864'  :'   m(yz)','201326593' :' 2/m(yz)','2097152'   :'  2(0+-)','1048576'   :'  m(0+-)',
1950    '3145729'   :'2/m(0+-)','8388608'   :'   2(xz)','4194304'   :'   m(xz)','12582913'  :' 2/m(xz)',
1951    '524288'    :'  2(+0-)','262144'    :'  m(+0-)','796433'    :'2/m(+0-)','1024'      :'   2(xy)',
1952    '512'       :'   m(xy)','1537'      :' 2/m(xy)','256'       :'  2(+-0)','128'       :'  m(+-0)',
1953    '385'       :'2/m(+-0)','76'        :'  mm2(x)','52'        :'  mm2(y)','42'        :'  mm2(z)',
1954    '135266336' :' mm2(yz)','69206048'  :'mm2(0+-)','8650760'   :' mm2(xz)','4718600'   :'mm2(+0-)',
1955    '1156'      :' mm2(xy)','772'       :'mm2(+-0)','82'        :'  222   ','136314944' :'  222(x)',
1956    '8912912'   :'  222(y)','1282'      :'  222(z)','127'       :'  mmm   ','204472417' :'  mmm(x)',
1957    '13369369'  :'  mmm(y)','1927'      :'  mmm(z)','33554496'  :'  4(100)','16777280'  :' -4(100)',
1958    '50331745'  :'4/m(100)','169869394' :'422(100)','84934738'  :'-42m 100','101711948' :'4mm(100)',
1959    '254804095' :'4/mmm100','536870928 ':'  4(010)','268435472' :' -4(010)','805306393' :'4/m (10)',
1960    '545783890' :'422(010)','272891986' :'-42m 010','541327412' :'4mm(010)','818675839' :'4/mmm010',
1961    '2050'      :'  4(001)','4098'      :' -4(001)','6151'      :'4/m(001)','3410'      :'422(001)',
1962    '4818'      :'-42m 001','2730'      :'4mm(001)','8191'      :'4/mmm001','8192'      :'  3(111)',
1963    '8193'      :' -3(111)','2629888'   :' 32(111)','1319040'   :' 3m(111)','3940737'   :'-3m(111)',
1964    '32768'     :'  3(+--)','32769'     :' -3(+--)','10519552'  :' 32(+--)','5276160'   :' 3m(+--)',
1965    '15762945'  :'-3m(+--)','65536'     :'  3(-+-)','65537'     :' -3(-+-)','134808576' :' 32(-+-)',
1966    '67437056'  :' 3m(-+-)','202180097' :'-3m(-+-)','131072'    :'  3(--+)','131073'    :' -3(--+)',
1967    '142737664' :' 32(--+)','71434368'  :' 3m(--+)','214040961' :'-3m(--+)','237650'    :'   23   ',
1968    '237695'    :'   m3   ','715894098' :'   432  ','358068946' :'  -43m  ','1073725439':'   m3m  ',
1969    '68157504'  :' mm2d100','4456464'   :' mm2d010','642'       :' mm2d001','153092172' :'-4m2 100',
1970    '277348404' :'-4m2 010','5418'      :'-4m2 001','1075726335':'  6/mmm ','1074414420':'-6m2 100',
1971    '1075070124':'-6m2 120','1075069650':'   6mm  ','1074414890':'   622  ','1073758215':'   6/m  ',
1972    '1073758212':'   -6   ','1073758210':'    6   ','1073759865':'-3m(100)','1075724673':'-3m(120)',
1973    '1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
1974    '1073758209':'   -3   ','1073758208':'    3   ','1074135143':'mmm(100)','1075314719':'mmm(010)',
1975    '1073743751':'mmm(110)','1074004034':' mm2z100','1074790418':' mm2z010','1073742466':' mm2z110',
1976    '1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
1977    '1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
1978    '1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
1979    '1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
1980    '1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
1981    '1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
1982    '1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
1983    '1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':'   -1   ',
1984    }
1985
1986NXUPQsym = {
1987    '    1   ':(28,29,28,28),'   -1   ':( 1,29,28, 0),'    2(x)':(12,18,12,25),'    m(x)':(25,18,12,25),
1988    '  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),
1989    '    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),
1990    '   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),
1991    '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),
1992    '  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),
1993    '   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),
1994    '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),
1995    ' mm2(yz)':(10,13, 0,-1),'mm2(0+-)':(11,13, 0,-1),' mm2(xz)':( 8,12, 0,-1),'mm2(+0-)':( 9,12, 0,-1),
1996    ' mm2(xy)':( 6,11, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'  222   ':( 1,10, 0,-1),'  222(x)':( 1,13, 0,-1),
1997    '  222(y)':( 1,12, 0,-1),'  222(z)':( 1,11, 0,-1),'  mmm   ':( 1,10, 0,-1),'  mmm(x)':( 1,13, 0,-1),
1998    '  mmm(y)':( 1,12, 0,-1),'  mmm(z)':( 1,11, 0,-1),'  4(100)':(12, 4,12, 0),' -4(100)':( 1, 4,12, 0),
1999    '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),
2000    '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),
2001    '422(010)':( 1, 3, 0,-1),'-42m 010':( 1, 3, 0,-1),'4mm(010)':(13, 3, 0,-1),'4/mmm010':(1, 3, 0,-1,),
2002    '  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),
2003    '-42m 001':( 1, 2, 0,-1),'4mm(001)':(14, 2, 0,-1),'4/mmm001':( 1, 2, 0,-1),'  3(111)':( 2, 5, 2, 0),
2004    ' -3(111)':( 1, 5, 2, 0),' 32(111)':( 1, 5, 0, 2),' 3m(111)':( 2, 5, 0, 2),'-3m(111)':( 1, 5, 0,-1),
2005    '  3(+--)':( 5, 8, 5, 0),' -3(+--)':( 1, 8, 5, 0),' 32(+--)':( 1, 8, 0, 5),' 3m(+--)':( 5, 8, 0, 5),
2006    '-3m(+--)':( 1, 8, 0,-1),'  3(-+-)':( 4, 7, 4, 0),' -3(-+-)':( 1, 7, 4, 0),' 32(-+-)':( 1, 7, 0, 4),
2007    ' 3m(-+-)':( 4, 7, 0, 4),'-3m(-+-)':( 1, 7, 0,-1),'  3(--+)':( 3, 6, 3, 0),' -3(--+)':( 1, 6, 3, 0),
2008    ' 32(--+)':( 1, 6, 0, 3),' 3m(--+)':( 3, 6, 0, 3),'-3m(--+)':( 1, 6, 0,-1),'   23   ':( 1, 1, 0, 0),
2009    '   m3   ':( 1, 1, 0, 0),'   432  ':( 1, 1, 0, 0),'  -43m  ':( 1, 1, 0, 0),'   m3m  ':( 1, 1, 0, 0),
2010    ' mm2d100':(12,13, 0,-1),' mm2d010':(13,12, 0,-1),' mm2d001':(14,11, 0,-1),'-4m2 100':( 1, 4, 0,-1),
2011    '-4m2 010':( 1, 3, 0,-1),'-4m2 001':( 1, 2, 0,-1),'  6/mmm ':( 1, 9, 0,-1),'-6m2 100':( 1, 9, 0,-1),
2012    '-6m2 120':( 1, 9, 0,-1),'   6mm  ':(14, 9, 0,-1),'   622  ':( 1, 9, 0,-1),'   6/m  ':( 1, 9,14,-1),
2013    '   -6   ':( 1, 9,14, 0),'    6   ':(14, 9,14, 0),'-3m(100)':( 1, 9, 0,-1),'-3m(120)':( 1, 9, 0,-1),
2014    ' 3m(100)':(14, 9, 0,14),' 3m(120)':(14, 9, 0,14),' 32(100)':( 1, 9, 0,14),' 32(120)':( 1, 9, 0,14),
2015    '   -3   ':( 1, 9,14, 0),'    3   ':(14, 9,14, 0),'mmm(100)':( 1,14, 0,-1),'mmm(010)':( 1,15, 0,-1),
2016    'mmm(110)':( 1,11, 0,-1),' mm2z100':(14,14, 0,-1),' mm2z010':(14,15, 0,-1),' mm2z110':(14,11, 0,-1),
2017    'mm2(100)':(12,14, 0,-1),'mm2(010)':(13,15, 0,-1),'mm2(110)':( 6,11, 0,-1),'mm2(120)':(15,14, 0,-1),
2018    'mm2(210)':(16,15, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'222(100)':( 1,14, 0,-1),'222(010)':( 1,15, 0,-1),
2019    '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),
2020    '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),
2021    ' m(001) ':(23,16,14,23),' m(100) ':(26,25,12,26),' m(010) ':(27,28,13,27),' m(110) ':(18,19, 6,18),
2022    ' m(120) ':(24,27,15,24),' m(210) ':(25,26,16,25),' m(+-0) ':(17,20, 7,17),' 2(001) ':(14,16,14,23),
2023    ' 2(100) ':(12,25,12,26),' 2(010) ':(13,28,13,27),' 2(110) ':( 6,19, 6,18),' 2(120) ':(15,27,15,24),
2024    ' 2(210) ':(16,26,16,25),' 2(+-0) ':( 7,20, 7,17),'   -1   ':( 1,29,28, 0)
2025    }
2026       
2027CSxinel = [[],      # 0th empty - indices are Fortran style
2028    [[0,0,0],[ 0.0, 0.0, 0.0]],      #1  0  0  0
2029    [[1,1,1],[ 1.0, 1.0, 1.0]],      #2  X  X  X
2030    [[1,1,1],[ 1.0, 1.0,-1.0]],      #3  X  X -X
2031    [[1,1,1],[ 1.0,-1.0, 1.0]],      #4  X -X  X
2032    [[1,1,1],[ 1.0,-1.0,-1.0]],      #5 -X  X  X
2033    [[1,1,0],[ 1.0, 1.0, 0.0]],      #6  X  X  0
2034    [[1,1,0],[ 1.0,-1.0, 0.0]],      #7  X -X  0
2035    [[1,0,1],[ 1.0, 0.0, 1.0]],      #8  X  0  X
2036    [[1,0,1],[ 1.0, 0.0,-1.0]],      #9  X  0 -X
2037    [[0,1,1],[ 0.0, 1.0, 1.0]],      #10  0  Y  Y
2038    [[0,1,1],[ 0.0, 1.0,-1.0]],      #11 0  Y -Y
2039    [[1,0,0],[ 1.0, 0.0, 0.0]],      #12  X  0  0
2040    [[0,1,0],[ 0.0, 1.0, 0.0]],      #13  0  Y  0
2041    [[0,0,1],[ 0.0, 0.0, 1.0]],      #14  0  0  Z
2042    [[1,1,0],[ 1.0, 2.0, 0.0]],      #15  X 2X  0
2043    [[1,1,0],[ 2.0, 1.0, 0.0]],      #16 2X  X  0
2044    [[1,1,2],[ 1.0, 1.0, 1.0]],      #17  X  X  Z
2045    [[1,1,2],[ 1.0,-1.0, 1.0]],      #18  X -X  Z
2046    [[1,2,1],[ 1.0, 1.0, 1.0]],      #19  X  Y  X
2047    [[1,2,1],[ 1.0, 1.0,-1.0]],      #20  X  Y -X
2048    [[1,2,2],[ 1.0, 1.0, 1.0]],      #21  X  Y  Y
2049    [[1,2,2],[ 1.0, 1.0,-1.0]],      #22  X  Y -Y
2050    [[1,2,0],[ 1.0, 1.0, 0.0]],      #23  X  Y  0
2051    [[1,0,2],[ 1.0, 0.0, 1.0]],      #24  X  0  Z
2052    [[0,1,2],[ 0.0, 1.0, 1.0]],      #25  0  Y  Z
2053    [[1,1,2],[ 1.0, 2.0, 1.0]],      #26  X 2X  Z
2054    [[1,1,2],[ 2.0, 1.0, 1.0]],      #27 2X  X  Z
2055    [[1,2,3],[ 1.0, 1.0, 1.0]],      #28  X  Y  Z
2056    ]
2057
2058CSuinel = [[],      # 0th empty - indices are Fortran style
2059    [[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
2060    [[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
2061    [[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
2062    [[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
2063    [[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
2064    [[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
2065    [[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
2066    [[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
2067    [[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
2068    [[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
2069    [[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
2070    [[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
2071    [[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
2072    [[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
2073    [[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
2074    [[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
2075    [[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
2076    [[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
2077    [[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
2078    [[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
2079    [[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
2080    [[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
2081    [[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
2082    [[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
2083    [[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
2084    [[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
2085    [[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
2086    [[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
2087    [[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
2088    ]
2089   
2090################################################################################
2091#### Site symmetry routines
2092################################################################################
2093   
2094def GetOprPtrName(key):
2095    'Needs a doc string'
2096    return OprPtrName[key]
2097
2098def GetOprName(key):
2099    'Needs a doc string'
2100    return OprName[key]
2101
2102def GetKNsym(key):
2103    'Needs a doc string'
2104    try:
2105        return KNsym[key]
2106    except KeyError:
2107        return 'sp'
2108
2109def GetNXUPQsym(siteSym):
2110    '''       
2111    The codes XUPQ are for lookup of symmetry constraints for position(X), thermal parm(U) & magnetic moments (P & Q)
2112    '''
2113    return NXUPQsym[siteSym]
2114
2115def GetCSxinel(siteSym): 
2116    "returns Xyz terms, multipliers, GUI flags"
2117    indx = GetNXUPQsym(siteSym)
2118    return CSxinel[indx[0]]
2119   
2120def GetCSuinel(siteSym):
2121    "returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
2122    indx = GetNXUPQsym(siteSym)
2123    return CSuinel[indx[1]]
2124   
2125def GetCSpqinel(siteSym,SpnFlp,dupDir): 
2126    "returns Mxyz terms, multipliers, GUI flags"
2127    CSI = [[1,2,3],[1.0,1.0,1.0]]
2128    for opr in dupDir:
2129        if '-1' in siteSym and SpnFlp[len(SpnFlp)//2] < 0:
2130            return [[0,0,0],[0.,0.,0.]]
2131        indx = GetNXUPQsym(opr)
2132        if SpnFlp[dupDir[opr]] > 0.:
2133            csi = CSxinel[indx[2]]  #P
2134        else:
2135            csi = CSxinel[indx[3]]  #Q
2136        if not len(csi):
2137            return [[0,0,0],[0.,0.,0.]]
2138        for kcs in [0,1,2]:
2139            if csi[0][kcs] == 0 and CSI[0][kcs] != 0:
2140                jcs = CSI[0][kcs]
2141                for ics in [0,1,2]:
2142                    if CSI[0][ics] == jcs:
2143                        CSI[0][ics] = 0
2144                        CSI[1][ics] = 0.
2145                    elif CSI[0][ics] > jcs:
2146                        CSI[0][ics] = CSI[0][jcs]-1
2147            elif CSI[0][kcs] == csi[0][kcs] and CSI[1][kcs] != csi[1][kcs]:
2148                CSI[1][kcs] = csi[1][kcs]
2149            elif CSI[0][kcs] > csi[0][kcs]:
2150                CSI[0][kcs] = min(CSI[0][kcs],csi[0][kcs])
2151                if CSI[1][kcs] == 1.:
2152                    CSI[1][kcs] = csi[1][kcs]
2153    return CSI
2154   
2155def getTauT(tau,sop,ssop,XYZ):
2156    ssopinv = nl.inv(ssop[0])
2157    mst = ssopinv[3][:3]
2158    epsinv = ssopinv[3][3]
2159    sdet = nl.det(sop[0])
2160    ssdet = nl.det(ssop[0])
2161    dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
2162    dT = 1.0
2163    if np.any(dtau%.5):
2164        dT = np.tan(np.pi*np.sum(dtau%.5))
2165    tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3])
2166    return sdet,ssdet,dtau,dT,tauT
2167   
2168def OpsfromStringOps(A,SGData,SSGData):
2169    SGOps = SGData['SGOps']
2170    SSGOps = SSGData['SSGOps']
2171    Ax = A.split('+')
2172    Ax[0] = int(Ax[0])
2173    iC = 1
2174    if Ax[0] < 0:
2175        iC = -1
2176    Ax[0] = abs(Ax[0])
2177    nA = Ax[0]%100-1
2178    return SGOps[nA],SSGOps[nA],iC
2179   
2180def GetSSfxuinel(waveType,nH,XYZ,SGData,SSGData,debug=False):
2181   
2182    def orderParms(CSI):
2183        parms = [0,]
2184        for csi in CSI:
2185            for i in [0,1,2]:
2186                if csi[i] not in parms:
2187                    parms.append(csi[i])
2188        for csi in CSI:
2189            for i in [0,1,2]:
2190                csi[i] = parms.index(csi[i])
2191        return CSI
2192       
2193    def fracCrenel(tau,Toff,Twid):
2194        Tau = (tau-Toff[:,np.newaxis])%1.
2195        A = np.where(Tau<Twid[:,np.newaxis],1.,0.)
2196        return A
2197       
2198    def fracFourier(tau,nH,fsin,fcos):
2199        SA = np.sin(2.*nH*np.pi*tau)
2200        CB = np.cos(2.*nH*np.pi*tau)
2201        A = SA[np.newaxis,np.newaxis,:]*fsin[:,:,np.newaxis]
2202        B = CB[np.newaxis,np.newaxis,:]*fcos[:,:,np.newaxis]
2203        return A+B
2204       
2205    def posFourier(tau,nH,psin,pcos):
2206        SA = np.sin(2*nH*np.pi*tau)
2207        CB = np.cos(2*nH*np.pi*tau)
2208        A = SA[np.newaxis,np.newaxis,:]*psin[:,:,np.newaxis]
2209        B = CB[np.newaxis,np.newaxis,:]*pcos[:,:,np.newaxis]
2210        return A+B   
2211
2212    def posSawtooth(tau,Toff,slopes):
2213        Tau = (tau-Toff)%1.
2214        A = slopes[:,np.newaxis]*Tau
2215        return A
2216   
2217    def posZigZag(tau,Tmm,XYZmax):
2218        DT = Tmm[1]-Tmm[0]
2219        slopeUp = 2.*XYZmax/DT
2220        slopeDn = 2.*XYZmax/(1.-DT)
2221        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])
2222        return A
2223
2224    def posBlock(tau,Tmm,XYZmax):
2225        A = np.array([np.where(Tmm[0] < t <= Tmm[1],XYZmax,-XYZmax) for t in tau])
2226        return A
2227       
2228    def DoFrac():
2229        delt2 = np.eye(2)*0.001
2230        FSC = np.ones(2,dtype='i')
2231        CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
2232        if 'Crenel' in waveType:
2233            dF = np.zeros_like(tau)
2234        else:
2235            dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2236        dFT = np.zeros_like(dF)
2237        dFTP = []
2238        for i in SdIndx:
2239            sop = Sop[i]
2240            ssop = SSop[i]           
2241            sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2242            fsc = np.ones(2,dtype='i')
2243            if 'Crenel' in waveType:
2244                dFT = np.zeros_like(tau)
2245                fsc = [1,1]
2246            else:   #Fourier
2247                dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
2248                dFT = nl.det(sop[0])*dFT
2249                dFT = dFT[:,np.argsort(tauT)]
2250                dFT[0] *= ssdet
2251                dFT[1] *= sdet
2252                dFTP.append(dFT)
2253           
2254                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2255                    fsc = [1,1]
2256                    CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
2257                    FSC = np.zeros(2,dtype='i')
2258                    return CSI,dF,dFTP
2259                else:
2260                    for i in range(2):
2261                        if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
2262                            fsc[i] = 1
2263                        else:
2264                            fsc[i] = 0
2265                    FSC &= fsc
2266                    if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc)
2267        n = -1
2268        for i,F in enumerate(FSC):
2269            if F:
2270                n += 1
2271                CSI[0][i] = n+1
2272                CSI[1][i] = 1.0
2273       
2274        return CSI,dF,dFTP
2275       
2276    def DoXYZ():
2277        delt4 = np.ones(4)*0.001
2278        delt5 = np.ones(5)*0.001
2279        delt6 = np.eye(6)*0.001
2280        if 'Fourier' in waveType:
2281            dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2282              #3x6x12 modulated position array (X,Spos,tau)& force positive
2283            CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2284        elif waveType == 'Sawtooth':
2285            dX = posSawtooth(tau,delt4[0],delt4[1:])
2286            CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0]]),
2287                np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2288        elif waveType in ['ZigZag','Block']:
2289            if waveType == 'ZigZag':
2290                dX = posZigZag(tau,delt5[:2],delt5[2:])
2291            else:
2292                dX = posBlock(tau,delt5[:2],delt5[2:])
2293            CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2294                np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2295        XSC = np.ones(6,dtype='i')
2296        dXTP = []
2297        for i in SdIndx:
2298            sop = Sop[i]
2299            ssop = SSop[i]
2300            sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2301            xsc = np.ones(6,dtype='i')
2302            if 'Fourier' in waveType:
2303                dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2304            elif waveType == 'Sawtooth':
2305                dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2306            elif waveType == 'ZigZag':
2307                dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2308            elif waveType == 'Block':
2309                dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2310            dXT = np.inner(sop[0],dXT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2311            dXT = np.swapaxes(dXT,1,2)      # back to array(3x6x49)
2312            dXT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2313            dXTP.append(dXT)
2314            if waveType == 'Fourier':
2315                for i in range(3):
2316                    if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2317                        xsc[i] = 0
2318                    if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2319                        xsc[i+3] = 0
2320                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2321                    xsc[3:6] = 0
2322                    CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2323                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]                   
2324                    if '(x)' in siteSym:
2325                        CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2326                        if 'm' in siteSym and len(SdIndx) == 1:
2327                            CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2328                    elif '(y)' in siteSym:
2329                        CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2330                        if 'm' in siteSym and len(SdIndx) == 1:
2331                            CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2332                    elif '(z)' in siteSym:
2333                        CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2334                        if 'm' in siteSym and len(SdIndx) == 1:
2335                            CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2336                if '4/mmm' in laue:
2337                    if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2338                        if '(xy)' in siteSym:
2339                            CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2340                            CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2341                    if '(xy)' in siteSym or '(+-0)' in siteSym:
2342                        mul = 1
2343                        if '(+-0)' in siteSym:
2344                            mul = -1
2345                        if np.allclose(dX[0,0,:],dXT[1,0,:]):
2346                            CSI[0][3:5] = [[11,0,0],[11,0,0]]
2347                            CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2348                            xsc[3:5] = 0
2349                        if np.allclose(dX[0,3,:],dXT[0,4,:]):
2350                            CSI[0][:2] = [[12,0,0],[12,0,0]]
2351                            CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2352                            xsc[:2] = 0
2353            XSC &= xsc
2354            if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2355        if waveType == 'Fourier':
2356            n = -1
2357            if debug: print (XSC)
2358            for i,X in enumerate(XSC):
2359                if X:
2360                    n += 1
2361                    CSI[0][i][0] = n+1
2362                    CSI[1][i][0] = 1.0
2363       
2364        return CSI,dX,dXTP
2365       
2366    def DoUij():
2367        tau = np.linspace(0,1,49,True)
2368        delt12 = np.eye(12)*0.0001
2369        dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
2370        CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2371        USC = np.ones(12,dtype='i')
2372        dUTP = []
2373        for i in SdIndx:
2374            sop = Sop[i]
2375            ssop = SSop[i]
2376            sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2377            usc = np.ones(12,dtype='i')
2378            dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
2379            dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3)    #convert dUT to 12x49x3x3
2380            dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2381            dUT = np.array(U2Uij(dUijT))    #convert to 6x12x49
2382            dUT = dUT[:,:,np.argsort(tauT)]
2383            dUT[:,:6,:] *=(ssdet*sdet)
2384            dUTP.append(dUT)
2385            if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2386                CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2387                [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2388                [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2389                [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2390                if 'mm2(x)' in siteSym:
2391                    CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2392                    USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2393                elif '(xy)' in siteSym:
2394                    CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2395                        [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2396                    CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2397                    USC = [1,1,1,1,1,1,1,1,1,1,1,1]                             
2398                elif '(x)' in siteSym:
2399                    CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2400                elif '(y)' in siteSym:
2401                    CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2402                elif '(z)' in siteSym:
2403                    CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2404                for i in range(6):
2405                    if not USC[i]:
2406                        CSI[0][i] = [0,0,0]
2407                        CSI[1][i] = [0.,0.,0.]
2408                        CSI[0][i+6] = [0,0,0]
2409                        CSI[1][i+6] = [0.,0.,0.]
2410            else:                       
2411                for i in range(6):
2412                    if not np.allclose(dU[i,i,:],dUT[i,i,:]):  #sin part
2413                        usc[i] = 0
2414                    if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]):   #cos part
2415                        usc[i+6] = 0
2416                if np.any(dUT[1,0,:]):
2417                    if '4/m' in siteSym:
2418                        CSI[0][6:8] = [[12,0,0],[12,0,0]]
2419                        if ssop[1][3]:
2420                            CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2421                            usc[9] = 1
2422                        else:
2423                            CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2424                            usc[9] = 0
2425                    elif '4' in siteSym:
2426                        CSI[0][6:8] = [[12,0,0],[12,0,0]]
2427                        CSI[0][:2] = [[11,0,0],[11,0,0]]
2428                        if ssop[1][3]:
2429                            CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2430                            CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2431                            usc[2] = 0
2432                            usc[8] = 0
2433                            usc[3] = 1
2434                            usc[9] = 1
2435                        else:
2436                            CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2437                            CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2438                            usc[2] = 1
2439                            usc[8] = 1
2440                            usc[3] = 0               
2441                            usc[9] = 0
2442                    elif 'xy' in siteSym or '+-0' in siteSym:
2443                        if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2444                            CSI[0][4:6] = [[12,0,0],[12,0,0]]
2445                            CSI[0][6:8] = [[11,0,0],[11,0,0]]
2446                            CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2447                            CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2448                            usc[4:6] = 0
2449                            usc[6:8] = 0
2450                       
2451                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2452            USC &= usc
2453        if debug: print (USC)
2454        if not np.any(dtau%.5):
2455            n = -1
2456            for i,U in enumerate(USC):
2457                if U:
2458                    n += 1
2459                    CSI[0][i][0] = n+1
2460                    CSI[1][i][0] = 1.0
2461
2462        return CSI,dU,dUTP
2463       
2464    if debug: print ('super space group: '+SSGData['SSpGrp'])
2465    CSI = {'Sfrac':[[[1,0],[2,0]],[[1.,0.],[1.,0.]]],
2466        'Spos':[[[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],
2467            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]],    #sin & cos
2468        'Sadp':[[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2469            [7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],
2470            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2471            [1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]],
2472        'Smag':[[[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],
2473            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]],}
2474    xyz = np.array(XYZ)%1.
2475    SGOps = copy.deepcopy(SGData['SGOps'])
2476    laue = SGData['SGLaue']
2477    siteSym = SytSym(XYZ,SGData)[0].strip()
2478    if debug: print ('siteSym: '+siteSym)
2479    if siteSym == '1':   #"1" site symmetry
2480        if debug:
2481            return CSI,None,None,None,None
2482        else:
2483            return CSI
2484    elif siteSym == '-1':   #"-1" site symmetry
2485        CSI['Sfrac'][0] = [[1,0],[0,0]]
2486        CSI['Spos'][0] = [[1,0,0],[2,0,0],[3,0,0], [0,0,0],[0,0,0],[0,0,0]]
2487        CSI['Sadp'][0] = [[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0], 
2488        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]]
2489        if debug:
2490            return CSI,None,None,None,None
2491        else:
2492            return CSI
2493    SSGOps = copy.deepcopy(SSGData['SSGOps'])
2494    #expand ops to include inversions if any
2495    if SGData['SGInv']:
2496        for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
2497            SGOps.append([-op[0],-op[1]%1.])
2498            SSGOps.append([-sop[0],-sop[1]%1.])
2499    #build set of sym ops around special position       
2500    SSop = []
2501    Sop = []
2502    Sdtau = []
2503    for iop,Op in enumerate(SGOps):         
2504        nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
2505        if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
2506            SSop.append(SSGOps[iop])
2507            Sop.append(SGOps[iop])
2508            ssopinv = nl.inv(SSGOps[iop][0])
2509            mst = ssopinv[3][:3]
2510            epsinv = ssopinv[3][3]
2511            Sdtau.append(np.sum(mst*(XYZ-SGOps[iop][1])-epsinv*SSGOps[iop][1][3]))
2512    SdIndx = np.argsort(np.array(Sdtau))     # just to do in sensible order
2513    if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
2514    #setup displacement arrays
2515    tau = np.linspace(-1,1,49,True)
2516    #make modulation arrays - one parameter at a time
2517    #site fractions
2518    CSI['Sfrac'],dF,dFTP = DoFrac()
2519    #positions
2520    CSI['Spos'],dX,dXTP = DoXYZ()       
2521    #anisotropic thermal motion
2522    CSI['Sadp'],dU,dUTP = DoUij()
2523    CSI['Spos'][0] = orderParms(CSI['Spos'][0])
2524    CSI['Sadp'][0] = orderParms(CSI['Sadp'][0])           
2525    if debug:
2526        return CSI,tau,[dF,dFTP],[dX,dXTP],[dU,dUTP]
2527    else:
2528        return CSI
2529   
2530def MustrainNames(SGData):
2531    'Needs a doc string'
2532    laue = SGData['SGLaue']
2533    uniq = SGData['SGUniq']
2534    if laue in ['m3','m3m']:
2535        return ['S400','S220']
2536    elif laue in ['6/m','6/mmm','3m1']:
2537        return ['S400','S004','S202']
2538    elif laue in ['31m','3']:
2539        return ['S400','S004','S202','S211']
2540    elif laue in ['3R','3mR']:
2541        return ['S400','S220','S310','S211']
2542    elif laue in ['4/m','4/mmm']:
2543        return ['S400','S004','S220','S022']
2544    elif laue in ['mmm']:
2545        return ['S400','S040','S004','S220','S202','S022']
2546    elif laue in ['2/m']:
2547        SHKL = ['S400','S040','S004','S220','S202','S022']
2548        if uniq == 'a':
2549            SHKL += ['S013','S031','S211']
2550        elif uniq == 'b':
2551            SHKL += ['S301','S103','S121']
2552        elif uniq == 'c':
2553            SHKL += ['S130','S310','S112']
2554        return SHKL
2555    else:
2556        SHKL = ['S400','S040','S004','S220','S202','S022']
2557        SHKL += ['S310','S103','S031','S130','S301','S013']
2558        SHKL += ['S211','S121','S112']
2559        return SHKL
2560       
2561def HStrainVals(HSvals,SGData):
2562    laue = SGData['SGLaue']
2563    uniq = SGData['SGUniq']
2564    DIJ = np.zeros(6)
2565    if laue in ['m3','m3m']:
2566        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
2567    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2568        DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
2569    elif laue in ['3R','3mR']:
2570        DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
2571    elif laue in ['4/m','4/mmm']:
2572        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
2573    elif laue in ['mmm']:
2574        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2575    elif laue in ['2/m']:
2576        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2577        if uniq == 'a':
2578            DIJ[5] = HSvals[3]
2579        elif uniq == 'b':
2580            DIJ[4] = HSvals[3]
2581        elif uniq == 'c':
2582            DIJ[3] = HSvals[3]
2583    else:
2584        DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
2585    return DIJ
2586
2587def HStrainNames(SGData):
2588    'Needs a doc string'
2589    laue = SGData['SGLaue']
2590    uniq = SGData['SGUniq']
2591    if laue in ['m3','m3m']:
2592        return ['D11','eA']         #add cubic strain term
2593    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2594        return ['D11','D33']
2595    elif laue in ['3R','3mR']:
2596        return ['D11','D12']
2597    elif laue in ['4/m','4/mmm']:
2598        return ['D11','D33']
2599    elif laue in ['mmm']:
2600        return ['D11','D22','D33']
2601    elif laue in ['2/m']:
2602        Dij = ['D11','D22','D33']
2603        if uniq == 'a':
2604            Dij += ['D23']
2605        elif uniq == 'b':
2606            Dij += ['D13']
2607        elif uniq == 'c':
2608            Dij += ['D12']
2609        return Dij
2610    else:
2611        Dij = ['D11','D22','D33','D12','D13','D23']
2612        return Dij
2613   
2614def MustrainCoeff(HKL,SGData):
2615    'Needs a doc string'
2616    #NB: order of terms is the same as returned by MustrainNames
2617    laue = SGData['SGLaue']
2618    uniq = SGData['SGUniq']
2619    h,k,l = HKL
2620    Strm = []
2621    if laue in ['m3','m3m']:
2622        Strm.append(h**4+k**4+l**4)
2623        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
2624    elif laue in ['6/m','6/mmm','3m1']:
2625        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
2626        Strm.append(l**4)
2627        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
2628    elif laue in ['31m','3']:
2629        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
2630        Strm.append(l**4)
2631        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
2632        Strm.append(4.0*h*k*l*(h+k))
2633    elif laue in ['3R','3mR']:
2634        Strm.append(h**4+k**4+l**4)
2635        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
2636        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
2637        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
2638    elif laue in ['4/m','4/mmm']:
2639        Strm.append(h**4+k**4)
2640        Strm.append(l**4)
2641        Strm.append(3.0*(h*k)**2)
2642        Strm.append(3.0*((h*l)**2+(k*l)**2))
2643    elif laue in ['mmm']:
2644        Strm.append(h**4)
2645        Strm.append(k**4)
2646        Strm.append(l**4)
2647        Strm.append(3.0*(h*k)**2)
2648        Strm.append(3.0*(h*l)**2)
2649        Strm.append(3.0*(k*l)**2)
2650    elif laue in ['2/m']:
2651        Strm.append(h**4)
2652        Strm.append(k**4)
2653        Strm.append(l**4)
2654        Strm.append(3.0*(h*k)**2)
2655        Strm.append(3.0*(h*l)**2)
2656        Strm.append(3.0*(k*l)**2)
2657        if uniq == 'a':
2658            Strm.append(2.0*k*l**3)
2659            Strm.append(2.0*l*k**3)
2660            Strm.append(4.0*k*l*h**2)
2661        elif uniq == 'b':
2662            Strm.append(2.0*l*h**3)
2663            Strm.append(2.0*h*l**3)
2664            Strm.append(4.0*h*l*k**2)
2665        elif uniq == 'c':
2666            Strm.append(2.0*h*k**3)
2667            Strm.append(2.0*k*h**3)
2668            Strm.append(4.0*h*k*l**2)
2669    else:
2670        Strm.append(h**4)
2671        Strm.append(k**4)
2672        Strm.append(l**4)
2673        Strm.append(3.0*(h*k)**2)
2674        Strm.append(3.0*(h*l)**2)
2675        Strm.append(3.0*(k*l)**2)
2676        Strm.append(2.0*k*h**3)
2677        Strm.append(2.0*h*l**3)
2678        Strm.append(2.0*l*k**3)
2679        Strm.append(2.0*h*k**3)
2680        Strm.append(2.0*l*h**3)
2681        Strm.append(2.0*k*l**3)
2682        Strm.append(4.0*k*l*h**2)
2683        Strm.append(4.0*h*l*k**2)
2684        Strm.append(4.0*k*h*l**2)
2685    return Strm
2686   
2687def Muiso2Shkl(muiso,SGData,cell):
2688    "this is to convert isotropic mustrain to generalized Shkls"
2689    import GSASIIlattice as G2lat
2690    A = G2lat.cell2AB(cell)[0]
2691   
2692    def minMus(Shkl,muiso,H,SGData,A):
2693        U = np.inner(A.T,H)
2694        S = np.array(MustrainCoeff(U,SGData))
2695        Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,np.newaxis]),axis=0))
2696        rad = np.sqrt(np.sum((Sum[:,np.newaxis]*H)**2,axis=1))
2697        return (muiso-rad)**2
2698       
2699    laue = SGData['SGLaue']
2700    PHI = np.linspace(0.,360.,60,True)
2701    PSI = np.linspace(0.,180.,60,True)
2702    X = np.outer(npsind(PHI),npsind(PSI))
2703    Y = np.outer(npcosd(PHI),npsind(PSI))
2704    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
2705    HKL = np.dstack((X,Y,Z))
2706    if laue in ['m3','m3m']:
2707        S0 = [1000.,1000.]
2708    elif laue in ['6/m','6/mmm','3m1']:
2709        S0 = [1000.,1000.,1000.]
2710    elif laue in ['31m','3']:
2711        S0 = [1000.,1000.,1000.,1000.]
2712    elif laue in ['3R','3mR']:
2713        S0 = [1000.,1000.,1000.,1000.]
2714    elif laue in ['4/m','4/mmm']:
2715        S0 = [1000.,1000.,1000.,1000.]
2716    elif laue in ['mmm']:
2717        S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
2718    elif laue in ['2/m']:
2719        S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
2720    else:
2721        S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000., 
2722            1000.,1000.,0.,0.,0.]
2723    S0 = np.array(S0)
2724    HKL = np.reshape(HKL,(-1,3))
2725    result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
2726    return result[0]
2727       
2728def PackRot(SGOps):
2729    IRT = []
2730    for ops in SGOps:
2731        M = ops[0]
2732        irt = 0
2733        for j in range(2,-1,-1):
2734            for k in range(2,-1,-1):
2735                irt *= 3
2736                irt += M[k][j]
2737        IRT.append(int(irt))
2738    return IRT
2739       
2740def SytSym(XYZ,SGData):
2741    '''
2742    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
2743
2744    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
2745    :param SGData: from SpcGroup
2746    :Returns: a two element tuple:
2747
2748     * The 1st element is a code for the site symmetry (see GetKNsym)
2749     * The 2nd element is the site multiplicity
2750
2751    '''
2752    Mult = 1
2753    Isym = 0
2754    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
2755        Isym = 1073741824
2756    Jdup = 0
2757    Ndup = 0
2758    dupDir = {}
2759    Xeqv = list(GenAtom(XYZ,SGData,True))
2760    IRT = PackRot(SGData['SGOps'])
2761    L = -1
2762    for ic,cen in enumerate(SGData['SGCen']):
2763        for invers in range(int(SGData['SGInv']+1)):
2764            for io,ops in enumerate(SGData['SGOps']):
2765                irtx = (1-2*invers)*IRT[io]
2766                L += 1
2767                if not Xeqv[L][1]:
2768                    Ndup = io
2769                    Jdup += 1
2770                    jx = GetOprPtrName(str(irtx))   #[KN table no,op name,KNsym ptr]
2771                    if jx[2] < 39:
2772                        px = GetOprName(str(irtx))
2773                        if px != '6643':    #skip Iden
2774                            dupDir[px] = io
2775                        Isym += 2**(jx[2]-1)
2776    if Isym == 1073741824: Isym = 0
2777    Mult = len(SGData['SGOps'])*len(SGData['SGCen'])*(int(SGData['SGInv'])+1)//Jdup
2778         
2779    return GetKNsym(str(Isym)),Mult,Ndup,dupDir
2780   
2781def ElemPosition(SGData):
2782    ''' Under development.
2783    Object here is to return a list of symmetry element types and locations suitable
2784    for say drawing them.
2785    So far I have the element type... getting all possible locations without lookup may be impossible!
2786    '''
2787    Inv = SGData['SGInv']
2788    eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
2789    # get operators & expand if centrosymmetric
2790    Ops = SGData['SGOps']
2791    opM = np.array([op[0].T for op in Ops])
2792    opT = np.array([op[1] for op in Ops])
2793    if Inv:
2794        opM = np.concatenate((opM,-opM))
2795        opT = np.concatenate((opT,-opT))
2796    opMT = list(zip(opM,opT))
2797    for M,T in opMT[1:]:        #skip I
2798        Dt = int(nl.det(M))
2799        Tr = int(np.trace(M))
2800        Dt = -(Dt-1)//2
2801        Es = eleSym[Tr][Dt]
2802        if Dt:              #rotation-inversion
2803            I = np.eye(3)
2804            if Tr == 1:     #mirrors/glides
2805                if np.any(T):       #glide
2806                    M2 = np.inner(M,M)
2807                    MT = np.inner(M,T)+T
2808                    print ('glide',Es,MT)
2809                    print (M2)
2810                else:               #mirror
2811                    print ('mirror',Es,T)
2812                    print (I-M)
2813                X = [-1,-1,-1]
2814            elif Tr == -3:  # pure inversion
2815                X = np.inner(nl.inv(I-M),T)
2816                print ('inversion',Es,X)
2817            else:           #other rotation-inversion
2818                M2 = np.inner(M,M)
2819                MT = np.inner(M,T)+T
2820                print ('rot-inv',Es,MT)
2821                print (M2)
2822                X = [-1,-1,-1]
2823        else:               #rotations
2824            print ('rotation',Es)
2825            X = [-1,-1,-1]
2826        #SymElements.append([Es,X])
2827       
2828    return #SymElements
2829   
2830def ApplyStringOps(A,SGData,X,Uij=[]):
2831    'Needs a doc string'
2832    SGOps = SGData['SGOps']
2833    SGCen = SGData['SGCen']
2834    Ax = A.split('+')
2835    Ax[0] = int(Ax[0])
2836    iC = 0
2837    if Ax[0] < 0:
2838        iC = 1
2839    Ax[0] = abs(Ax[0])
2840    nA = Ax[0]%100-1
2841    cA = Ax[0]//100
2842    Cen = SGCen[cA]
2843    M,T = SGOps[nA]
2844    if len(Ax)>1:
2845        cellA = Ax[1].split(',')
2846        cellA = np.array([int(a) for a in cellA])
2847    else:
2848        cellA = np.zeros(3)
2849    newX = Cen+(1-2*iC)*(np.inner(M,X).T+T)+cellA
2850    if len(Uij):
2851        U = Uij2U(Uij)
2852        U = np.inner(M,np.inner(U,M).T)
2853        newUij = U2Uij(U)
2854        return [newX,newUij]
2855    else:
2856        return newX
2857       
2858def ApplyStringOpsMom(A,SGData,Mom):
2859    'Needs a doc string'
2860    SGOps = SGData['SGOps']
2861    Ax = A.split('+')
2862    Ax[0] = int(Ax[0])
2863    Ax[0] = abs(Ax[0])
2864    nA = Ax[0]%100-1
2865    M,T = SGOps[nA]
2866    if len(Ax)>1:
2867        cellA = Ax[1].split(',')
2868        cellA = np.array([int(a) for a in cellA])
2869    else:
2870        cellA = np.zeros(3)
2871    if SGData['SGGray']:
2872        newMom = -(np.inner(Mom,M).T)*nl.det(M)
2873    else:
2874        newMom = -(np.inner(Mom,M).T)*SGData['SpnFlp'][nA-1]*nl.det(M)
2875    return newMom
2876       
2877def StringOpsProd(A,B,SGData):
2878    """
2879    Find A*B where A & B are in strings '-' + '100*c+n' + '+ijk'
2880    where '-' indicates inversion, c(>0) is the cell centering operator,
2881    n is operator number from SgOps and ijk are unit cell translations (each may be <0).
2882    Should return resultant string - C. SGData - dictionary using entries:
2883
2884       *  'SGCen': cell centering vectors [0,0,0] at least
2885       *  'SGOps': symmetry operations as [M,T] so that M*x+T = x'
2886
2887    """
2888    SGOps = SGData['SGOps']
2889    SGCen = SGData['SGCen']
2890    #1st split out the cell translation part & work on the operator parts
2891    Ax = A.split('+'); Bx = B.split('+')
2892    Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
2893    iC = 0
2894    if Ax[0]*Bx[0] < 0:
2895        iC = 1
2896    Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
2897    nA = Ax[0]%100-1;  nB = Bx[0]%100-1
2898    cA = Ax[0]//100;  cB = Bx[0]//100
2899    Cen = (SGCen[cA]+SGCen[cB])%1.0
2900    cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
2901    Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
2902    Mc = np.inner(Ma,Mb.T)
2903#    print Ma,Mb,Mc
2904    Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
2905#    print Ta,Tb,Tc
2906#    print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
2907    nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
2908    #now the cell translation part
2909    if len(Ax)>1:
2910        cellA = Ax[1].split(',')
2911        cellA = [int(a) for a in cellA]
2912    else:
2913        cellA = [0,0,0]
2914    if len(Bx)>1:
2915        cellB = Bx[1].split(',')
2916        cellB = [int(b) for b in cellB]
2917    else:
2918        cellB = [0,0,0]
2919    cellC = np.add(cellA,cellB)
2920    C = str(((nC+1)+(100*cC))*(1-2*iC))+'+'+ \
2921        str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
2922    return C
2923           
2924def U2Uij(U):
2925    #returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
2926    return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
2927   
2928def Uij2U(Uij):
2929    #returns the thermal motion tensor U from Uij as numpy array
2930    return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
2931
2932def StandardizeSpcName(spcgroup):
2933    '''Accept a spacegroup name where spaces may have not been used
2934    in the names according to the GSAS convention (spaces between symmetry
2935    for each axis) and return the space group name as used in GSAS
2936    '''
2937    rspc = spcgroup.replace(' ','').upper()
2938    # deal with rhombohedral and hexagonal setting designations
2939    rhomb = ''
2940    if rspc[-1:] == 'R':
2941        rspc = rspc[:-1]
2942        rhomb = ' R'
2943    gray = ''
2944    if "1'" in rspc:
2945        gray = " 1'"
2946        rspc = rspc.replace("1'",'')
2947    elif rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
2948        rspc = rspc[:-1]
2949    else:
2950        rspc = rspc.replace("'",'')
2951    # look for a match in the spacegroup lists
2952    for i in spglist.values():
2953        for spc in i:
2954            if rspc == spc.replace(' ','').upper():
2955                return spc+gray+rhomb
2956    # how about the post-2002 orthorhombic names?
2957    if rspc in sgequiv_2002_orthorhombic:
2958        return sgequiv_2002_orthorhombic[rspc]+gray
2959    else:
2960    # not found
2961        return ''
2962
2963spgbyNum = []
2964'''Space groups indexed by number'''
2965spgbyNum = [None,
2966        'P 1','P -1',                                                   #1-2
2967        'P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
2968        'C 2/m','P 2/c','P 21/c','C 2/c',                               #3-15
2969        'P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21',
2970        'C 2 2 21','C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
2971        'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
2972        'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
2973        'C m m 2','C m c 21','C c c 2',
2974        'A m m 2','A b m 2','A m a 2','A b a 2',
2975        'F m m 2','F d d 2','I m m 2','I b a 2','I m a 2',
2976        'P m m m','P n n n','P c c m','P b a n',
2977        'P m m a','P n n a','P m n a','P c c a','P b a m','P c c n',
2978        'P b c m','P n n m','P m m n','P b c n','P b c a','P n m a',
2979        'C m c m','C m c a','C m m m','C c c m','C m m a','C c c a',
2980        'F m m m', 'F d d d',
2981        'I m m m','I b a m','I b c a','I m m a',                        #16-74
2982        'P 4','P 41','P 42','P 43',
2983        'I 4','I 41',
2984        'P -4','I -4','P 4/m','P 42/m','P 4/n','P 42/n',
2985        'I 4/m','I 41/a',
2986        'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
2987        'P 42 21 2','P 43 2 2','P 43 21 2',
2988        'I 4 2 2','I 41 2 2',
2989        'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c','P 4 n c',
2990        'P 42 m c','P 42 b c',
2991        'I 4 m m','I 4 c m','I 41 m d','I 41 c d',
2992        'P -4 2 m','P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2',
2993        'P -4 c 2','P -4 b 2','P -4 n 2',
2994        'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d',
2995        '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',
2996        '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',
2997        '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',
2998        'P 42/n c m',
2999        'I 4/m m m','I 4/m c m','I 41/a m d','I 41/a c d',
3000        'P 3','P 31','P 32','R 3','P -3','R -3',
3001        'P 3 1 2','P 3 2 1','P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1',
3002        'R 3 2',
3003        'P 3 m 1','P 3 1 m','P 3 c 1','P 3 1 c',
3004        'R 3 m','R 3 c',
3005        'P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1',
3006        'R -3 m','R -3 c',                                               #75-167
3007        'P 6','P 61',
3008        'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3009        'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3010        'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3011        '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
3012        'P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m 3','P n 3',
3013        'F m -3','F d -3','I m -3',
3014        'P a -3','I a -3','P 4 3 2','P 42 3 2','F 4 3 2','F 41 3 2',
3015        'I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2','P -4 3 m',
3016        'F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
3017        'P m -3 m','P n -3 n','P m -3 n','P n -3 m',
3018        'F m -3 m','F m -3 c','F d -3 m','F d -3 c',
3019        'I m -3 m','I a -3 d',]                                       #195-230
3020spg2origins = {}
3021''' A dictionary of all spacegroups that have 2nd settings; the value is the
30221st --> 2nd setting transformation vector as X(2nd) = X(1st)-V, nonstandard ones are included.
3023'''
3024spg2origins = {
3025        'P n n n':[-.25,-.25,-.25],
3026        'P b a n':[-.25,-.25,0],'P n c b':[0,-.25,-.25],'P c n a':[-.25,0,-.25],
3027        'P m m n':[-.25,-.25,0],'P n m m':[0,-.25,-.25],'P m n m':[-.25,0,-.25],
3028        'C c c a':[0,-.25,-.25],'C c c b':[-.25,0,-.25],'A b a a':[-.25,0,-.25],
3029        'A c a a':[-.25,-.25,0],'B b c b':[-.25,-.25,0],'B b a b':[0,-.25,-.25],
3030        'F d d d':[-.125,-.125,-.125],
3031        'P 4/n':[-.25,-.25,0],'P 42/n':[-.25,-.25,-.25],'I 41/a':[0,-.25,-.125],
3032        '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],
3033        '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],
3034        'I 41/a m d':[0,.25,-.125],'I 41/a c d':[0,.25,-.125],
3035        'p n -3':[-.25,-.25,-.25],'F d -3':[-.125,-.125,-.125],'P n -3 n':[-.25,-.25,-.25],
3036        'P n -3 m':[-.25,-.25,-.25],'F d -3 m':[-.125,-.125,-.125],'F d -3 c':[-.375,-.375,-.375],
3037        'p n 3':[-.25,-.25,-.25],'F d 3':[-.125,-.125,-.125],'P n 3 n':[-.25,-.25,-.25],
3038        'P n 3 m':[-.25,-.25,-.25],'F d 3 m':[-.125,-.125,-.125],'F d - c':[-.375,-.375,-.375]}
3039spglist = {}
3040'''A dictionary of space groups as ordered and named in the pre-2002 International
3041Tables Volume A, except that spaces are used following the GSAS convention to
3042separate the different crystallographic directions.
3043Note that the symmetry codes here will recognize many non-standard space group
3044symbols with different settings. They are ordered by Laue group
3045'''
3046spglist = {
3047    'P1' : ('P 1','P -1',), # 1-2
3048    'C1' : ('C 1','C -1',),
3049    'P2/m': ('P 2','P 21','P m','P a','P c','P n',
3050        '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
3051    'C2/m':('C 2','C m','C c','C n',
3052        'C 2/m','C 2/c','C 2/n',),
3053    'A2/m':('A 2','A m','A a','A n',
3054        'A 2/m','A 2/a','A 2/n',),
3055   'Pmmm':('P 2 2 2',
3056        'P 2 2 21','P 21 2 2','P 2 21 2',
3057        'P 21 21 2','P 2 21 21','P 21 2 21',
3058        'P 21 21 21',
3059        'P m m 2','P 2 m m','P m 2 m',
3060        'P m c 21','P 21 m a','P b 21 m','P m 21 b','P c m 21','P 21 a m',
3061        'P c c 2','P 2 a a','P b 2 b',
3062        'P m a 2','P 2 m b','P c 2 m','P m 2 a','P b m 2','P 2 c m',
3063        'P c a 21','P 21 a b','P c 21 b','P b 21 a','P b c 21','P 21 c a',
3064        'P n c 2','P 2 n a','P b 2 n','P n 2 b','P c n 2','P 2 a n',
3065        'P m n 21','P 21 m n','P n 21 m','P m 21 n','P n m 21','P 21 n m',
3066        'P b a 2','P 2 c b','P c 2 a',
3067        'P n a 21','P 21 n b','P c 21 n','P n 21 a','P b n 21','P 21 c n',
3068        'P n n 2','P 2 n n','P n 2 n',
3069        'P m m m','P n n n',
3070        'P c c m','P m a a','P b m b',
3071        'P b a n','P n c b','P c n a',
3072        'P m m a','P b m m','P m c m','P m a m','P m m b','P c m m',
3073        'P n n a','P b n n','P n c n','P n a n','P n n b','P c n n',
3074        'P m n a','P b m n','P n c m','P m a n','P n m b','P c n m',
3075        'P c c a','P b a a','P b c b','P b a b','P c c b','P c a a',
3076        'P b a m','P m c b','P c m a',
3077        'P c c n','P n a a','P b n b',
3078        'P b c m','P m c a','P b m a','P c m b','P c a m','P m a b',
3079        'P n n m','P m n n','P n m n',
3080        'P m m n','P n m m','P m n m',
3081        'P b c n','P n c a','P b n a','P c n b','P c a n','P n a b',
3082        'P b c a','P c a b',
3083        'P n m a','P b n m','P m c n','P n a m','P m n b','P c m n',
3084        ),
3085    'Cmmm':('C 2 2 21','C 2 2 2','C m m 2',
3086        'C m c 21','C c m 21','C c c 2','C m 2 m','C 2 m m',
3087        'C m 2 a','C 2 m b','C c 2 m','C 2 c m','C c 2 a','C 2 c b',
3088        'C m c m','C c m m','C m c a','C c m b',
3089        'C m m m','C c c m','C m m a','C m m b','C c c a','C c c b',),
3090    'Ammm':('A 21 2 2','A 2 2 2','A 2 m m',
3091        'A 21 m a','A 21 a m','A 2 a a','A m 2 m','A m m 2',
3092        'A b m 2','A c 2 m','A m a 2','A m 2 a','A b a 2','A c 2 a',
3093        'A m m a','A m a m','A b m a','A c a m',
3094        'A m m m','A m a a','A b m m','A c m m','A c a a','A b a a',),
3095    'Bmmm':('B 2 21 2','B 2 2 2','B m 2 m',
3096        'B m 21 b','B b 21 m','B b 2 b','B m m 2','B 2 m m',
3097        'B 2 c m','B m a 2','B 2 m b','B b m 2','B 2 c b','B b a 2',
3098        'B b m m','B m m b','B b c m','B m a b',
3099        'B m m m','B b m b','B m a m','B m c m','B b a b','B b c b',),
3100    'Immm':('I 2 2 2','I 21 21 21',
3101        'I m m 2','I m 2 m','I 2 m m',
3102        'I b a 2','I 2 c b','I c 2 a',
3103        'I m a 2','I 2 m b','I c 2 m','I m 2 a','I b m 2','I 2 c m',
3104        'I m m m','I b a m','I m c b','I c m a',
3105        'I b c a','I c a b',
3106        'I m m a','I b m m ','I m c m','I m a m','I m m b','I c m m',),
3107    'Fmmm':('F 2 2 2','F m m m', 'F d d d',
3108        'F m m 2','F m 2 m','F 2 m m',
3109        'F d d 2','F d 2 d','F 2 d d',),
3110    'P4/mmm':('P 4','P 41','P 42','P 43','P -4','P 4/m','P 42/m','P 4/n','P 42/n',
3111        'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3112        'P 42 21 2','P 43 2 2','P 43 21 2','P 4 m m','P 4 b m','P 42 c m',
3113        'P 42 n m','P 4 c c','P 4 n c','P 42 m c','P 42 b c','P -4 2 m',
3114        'P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2',
3115        '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',
3116        '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',
3117        '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',
3118        'P 42/n c m',),
3119    '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',
3120        'I 4 c m','I 41 m d','I 41 c d',
3121        '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',
3122        'I 41/a m d','I 41/a c d'),
3123    'R3-H':('R 3','R -3','R 3 2','R 3 m','R 3 c','R -3 m','R -3 c',),
3124    'P6/mmm': ('P 3','P 31','P 32','P -3','P 3 1 2','P 3 2 1','P 31 1 2',
3125        'P 31 2 1','P 32 1 2','P 32 2 1', 'P 3 m 1','P 3 1 m','P 3 c 1',
3126        'P 3 1 c','P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1','P 6','P 61',
3127        'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3128        'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3129        'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3130        'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',),
3131    'Pm3m': ('P 2 3','P 21 3','P m 3','P n 3','P a 3','P 4 3 2','P 42 3 2',
3132        'P 43 3 2','P 41 3 2','P -4 3 m','P -4 3 n','P m 3 m','P n 3 n',
3133        'P m 3 n','P n 3 m','P n -3 n','P n -3 m','P m -3 n',),
3134    'Im3m':('I 2 3','I 21 3','I m -3','I a -3', 'I 4 3 2','I 41 3 2',
3135        'I -4 3 m', 'I -4 3 d','I m -3 m','I m 3 m','I a -3 d','I n -3 n'),
3136    'Fm3m':('F 2 3','F m -3','F d -3','F 4 3 2','F 41 3 2','F -4 3 m',
3137        'F -4 3 c','F m -3 m','F m 3 m','F m -3 c','F d -3 m','F d -3 c',),
3138}
3139sgequiv_2002_orthorhombic = {}
3140''' A dictionary of orthorhombic space groups that were renamed in the 2002 Volume A,
3141 along with the pre-2002 name. The e designates a double glide-plane
3142'''
3143sgequiv_2002_orthorhombic = {
3144        'AEM2':'A b m 2','B2EM':'B 2 c m','CM2E':'C m 2 a',
3145        'AE2M':'A c 2 m','BME2':'B m a 2','C2ME':'C 2 m b',
3146        'AEA2':'A b a 2','B2EB':'B 2 c b','CC2E':'C c 2 a',
3147        'AE2A':'A c 2 a','BBE2':'B b a 2','C2CE':'C 2 c b',
3148        'CMCE':'C m c a','AEMA':'A b m a','BBEM':'B b c m',
3149        'BMEB':'B m a b','CCME':'C c m b','AEAM':'A c a m',
3150        'CMME':'C m m a','AEMM':'A b m m','BMEM':'B m c m',
3151        'CCCE':'C c c a','AEAA':'A b a a','BBEB':'B b c b'}
3152ptssdict = {}
3153'''A dictionary of superspace group symbols allowed for each point group
3154(except cubics). Monoclinics are all b-unique setting.
3155'''
3156ptssdict = {
3157#1,2
3158    'P1':['(abg)',],'C1':['(abg)',],
3159#3-15
3160    'P2':['(a0g)','(a1/2g)','(0b0)','(0b0)s','(1/2b0)','(0b1/2)',],
3161    'C2':['(a0g)','(0b0)','(0b0)s','(0b1/2)',],
3162    'A2':['(a0g)','(0b0)','(0b0)s','(1/2b0)',],
3163    'Pm':['(a0g)','(a0g)s','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)',],
3164    'Cm':['(a0g)','(a0g)s','(0b0)','(0b1/2)',],
3165    'Am':['(a0g)','(a0g)s','(0b0)','(1/2b0)',],
3166    'P2/m':['(a0g)','(a0g)0s','(a1/2g)','(0b0)','(0b0)s0','(1/2b0)','(0b1/2)','(1/2b1/2)',],
3167    'C2/m':['(a0g)','(a0g)0s','(0b0)','(0b0)s0','(0b1/2)',],
3168    'A2/m':['(a0g)','(a0g)s0','(0b0)','(0b0)0s','(1/2b0)',],
3169#16-24
3170    'P222':['(00g)','(00g)00s','(01/2g)','(1/20g)','(1/21/2g)',
3171            '(a00)','(a00)s00','(a01/2)','(a1/20)','(a1/21/2)',
3172            '(0b0)','(0b0)0s0','(1/2b0)','(0b1/2)','(1/2b1/2)',],
3173    'C222':['(00g)','(00g)00s','(10g)','(10g)00s','(01g)','(01g)00s',
3174            '(a00)','(a00)s00','(a01/2)',
3175            '(0b0)','(0b0)0s0','(0b1/2)',],
3176    'A222':['(a00)','(a00)s00','(a10)','(a10)s00','(a01)','(a01)s00',
3177            '(0b0)','(0b0)0s0','(1/2b0)',
3178            '(00g)','(00g)00s','(1/20g)',],
3179    'B222':['(0b0)','(0b0)0s0','(1b0)','(1b0)0s0','(0b1)','(0b1)0s0',
3180            '(00g)','(00g)00s','(01/2g)',
3181            '(a00)','(a00)s00','(a1/20)',],
3182    'F222':['(00g)','(00g)00s','(10g)','(01g)',
3183            '(a00)','(a00)s00','(a10)','(a01)',
3184            '(0b0)','(0b0)0s0','(1b0)','(0b1)',],
3185    'I222':['(00g)','(00g)00s',
3186            '(a00)','(a00)s00',
3187            '(0b0)','(0b0)0s0',],
3188#25-46
3189    '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',
3190            '(a00)','(a00)0s0','(a1/20)','(a01/2)','(a01/2)0s0','(a1/21/2)','(1/21/2g)qq0',
3191            '(0b0)','(0b0)s00','(0b1/2)','(0b1/2)s00','(1/2b0)','(1/2b1/2)','(1/2b1/2)q00',],       
3192    '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',
3193            '(0b0)','(0b0)00s','(1/2b0)','(0b1/2)','(0b1/2)00s','(1/2b1/2)','(a1/21/2)0qq',
3194            '(00g)','(00g)0s0','(01/2g)','(01/2g)0s0','(1/20g)','(1/21/2g)','(1/21/2g)0q0',],
3195    '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',
3196            '(00g)','(00g)s00','(1/20g)','(01/2g)','(01/2g)s00','(1/21/2g)','(1/2b1/2)q0q',
3197            '(a00)','(a00)0s0','(a01/2)','(a01/2)0s0','(a1/20)','(a1/21/2)','(a1/21/2)00q',],
3198    'Cmm2':['(00g)','(00g)s0s','(00g)ss0','(10g)','(10g)s0s','(10g)ss0',
3199            '(a00)','(a00)0s0','(a01/2)','(a01/2)0s0',
3200            '(0b0)','(0b0)s00','(0b1/2)','(0b1/2)s00',],
3201    'C2mm':['(0b0)','(0b0)00s','(0b1/2)',],
3202    'Cm2m':['(0b0)','(0b0)0ss','(0b0)s0s','(0b0)ss0','(0b1/2)','(0b1/2)ss0',],
3203    'A2mm':['(a00)','(a00)ss0','(a00)0ss','(a10)','(a10)ss0','(a10)0ss',
3204            '(0b0)','(0b0)00s','(1/2b0)','(1/2b0)00s',
3205            '(00g)','(00g)0s0','(1/20g)','(1/20g)0s0',],
3206    'Am2m':['(00g)','(00g)s00','(1/20g)',],
3207    'Amm2':['(00g)','(00g)s0s','(00g)ss0','(00g)0ss','(1/20g)','(1/20g)0ss',],
3208    'Bm2m':['(0b0)','(0b0)0ss','(0b0)s0s','(0b1)','(0b1)0ss','(0b1)s0s',
3209            '(a00)','(a00)00s','(a1/20)','(a1/20)00s',
3210            '(00g)','(00g)s00','(01/2g)','(01/2g)s00',],
3211    'B2mm':['(a00)','(a00)0ss','(a00)s0s','(a00)ss0','(a1/20)','(a1/20)s0s',],
3212    'Bmm2':['(a00)','(a00)0s0','(a1/20)',],
3213    'Fmm2':['(00g)','(00g)s0s','(00g)0ss','(00g)ss0','(10g)','(10g)ss0','(10g)s0s','(01g)','(01g)ss0','(01g)0ss',
3214            '(a00)','(a00)0s0','(a01)','(a01)0s0',
3215            '(0b0)','(0b0)s00','(0b1)','(0b1)s00',],       
3216    'F2mm':['(a00)','(a00)ss0','(a00)s0s','(a00)0ss','(a10)','(a10)0ss','(a10)ss0','(a01)','(a01)0ss','(a01)s0s',
3217            '(0b0)','(0b0)00s','(1b0)','(1b0)00s',
3218            '(00g)','(00g)0s0','(10g)','(10g)0s0',],
3219    'Fm2m':['(0b0)','(0b0)0ss','(0b0)ss0','(0b0)s0s','(0b1)','(0b1)s0s','(0b1)0ss','(1b0)','(1b0)s0s','(1b0)ss0',
3220            '(00g)','(00g)s00','(01g)','(01g)s00',
3221            '(a00)','(a00)00s','(a10)','(a10)00s',],       
3222    'Imm2':['(00g)','(00g)ss0','(00g)s0s','(00g)0ss','(a00)','(a00)0s0','(0b0)','(0b0)s00',],
3223    'I2mm':['(a00)','(00g)0ss','(00g)ss0','(00g)s0s','(0b0)','(0b0)00s','(00g)','(00g)0s0',],
3224    'Im2m':['(0b0)','(0b0)s0s','(0b0)0ss','(0b0)ss0','(00g)','(00g)s00','(a00)','(a00)00s',],
3225#47-74
3226    'Pmmm':['(00g)','(00g)s00','(00g)0s0','(00g)ss0','(01/2g)','(01/2g)s00','(1/20g)','(1/20g)0s0','(1/21/2g)',
3227            '(a00)','(a00)0s0','(a00)00s','(a00)0ss','(a01/2)','(a01/2)0s0','(a1/20)','(a1/20)00s','(a1/21/2)',
3228            '(0b0)','(0b0)s00','(0b0)00s','(0b0)s0s','(1/2b0)','(1/2b0)00s','(0b1/2)','(0b1/2)s00','(1/2b1/2)',],
3229    'Cmmm':['(00g)','(00g)s00','(00g)ss0','(10g)','(10g)s00','(10g)ss0','(01g)','(01g)0s0','(01g)ss0',
3230            '(a00)','(a00)00s','(a00)0ss','(a00)0s0','(a01/2)','(a01/2)0s0',
3231            '(0b0)','(0b0)00s','(0b0)s0s','(0b0)s00','(0b1/2)','(0b1/2)s00',],
3232    'Ammm':['(a00)','(a00)0s0','(a00)0ss','(a10)','(a10)0s0','(a10)0ss','(a01)','(a01)0s0','(a01)0ss',
3233            '(0b0)','(0b0)00s','(0b0)s0s','(0b0)s00','(0b1/2)','(0b1/2)s00',
3234            '(00g)','(00g)s00','(00g)ss0','(00g)0s0','(1/20g)','(1/20g)0s0',],
3235    'Bmmm':['(0b0)','(0b0)00s','(0b0)s0s','(0b1)','(0b1)00s','(0b1)s0s','(1b0)','(1b0)00s','(1b0)s0s',
3236            '(a00)','(a00)0s0','(a00)0ss','(a00)00s','(a1/20)','(a1/20)00s',
3237            '(00g)','(00g)s00','(00g)ss0','(00g)0s0','(1/20g)','(1/20g)0s0',],
3238    'Fmmm':['(00g)','(00g)s00','(00g)ss0','(10g)','(10g)s00','(10g)0s0','(10g)ss0','(01g)','(01g)s00','(01g)0s0','(01g)ss0',
3239            '(a00)','(a00)0s0','(a00)0ss','(a10)','(a10)0s0','(a10)00s','(a10)0ss','(a01)','(a01)0s0','(a01)00s','(a01)0ss',
3240            '(0b0)','(0b0)s00','(0b0)s0s','(0b1)','(0b1)s00','(0b1)00s','(0b1)s0s','(1b0)','(1b0)s00','(1b0)00s','(1b0)s0s'],
3241#75-82
3242    'P4':['(00g)','(00g)q','(00g)s','(1/21/2g)','(1/21/2g)q',],
3243    'I4':['(00g)','(00g)q','(00g)s',],
3244    'P-4':['(00g)','(1/21/2g)',],
3245    'I-4':['(00g)',],
3246#83-89
3247    'P4/m':['(00g)','(00g)s0','(1/21/2g)',],
3248    'I4/m':['(00g)','(00g)s0',],
3249#90-98
3250    'P422':['(00g)','(00g)q00','(00g)s00','(1/21/2g)','(1/21/2g)q00',],
3251    'I422':['(00g)','(00g)q00','(00g)s00',],
3252#99-122
3253    'P4mm':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s','(1/21/2g)','(1/21/2g)0ss','(1/21/2g)qq0','(1/21/2g)qqs',],
3254    'I4mm':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',],
3255    'P-42m':['(00g)','(00g)0ss','(1/21/2g)','(1/21/2g)0ss',],
3256    'P-4m2':['(00g)','(00g)0s0','(1/21/2g)','(1/21/2g)0q0',],
3257    'I-4m2':['(00g)','(00g)0s0',],
3258    'I-42m':['(00g)','(00g)0ss',],
3259#123-142   
3260    'P4/mmm':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',
3261        '(1/21/2g)','(1/21/2g)s0s0','(1/21/2g)00ss','(1/21/2g)s00s',],
3262    'I4/mmm':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
3263#143-148
3264    'P 3':['(00g)','(00g)t','(1/31/3g)',],
3265    'R3':['(00g)','(00g)t',],
3266    'P-3':['(00g)','(1/31/3g)',],
3267    'R-3':['(00g)',],
3268#149-161
3269    'P312':['(00g)','(00g)t00','(1/31/3g)',],
3270    'P321':['(00g)','(00g)t00',],
3271    'R32':['(00g)','(00g)t0',],
3272    'P3m1':['(00g)','(00g)0s0',],
3273    'P31m':['(00g)','(00g)00s','(1/31/3g)','(1/31/3g)00s',],
3274    'R3m':['(00g)','(00g)0s',],
3275#162-167
3276    'P-31m':['(00g)','(00g)00s','(1/31/3g)','(1/31/3g)00s',],
3277    'P-3m1':['(00g)','(00g)0s0',],
3278    'R-3m':['(00g)','(00g)0s',],
3279#168-176
3280    'P6':['(00g)','(00g)h','(00g)t','(00g)s',],
3281    'P-6':['(00g)',],
3282    'P6/m':['(00g)','(00g)s0',],
3283#177-194
3284    'P622':['(00g)','(00g)h00','(00g)t00','(00g)s00',],
3285    'P6mm':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',],
3286    'P-6m2':['(00g)','(00g)0s0',],
3287    'P-62m':['(00g)','(00g)00s',],
3288    'P6/mmm':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
3289        }
3290
3291ssdict = {}
3292
3293#'A few non-standard space groups for test use'
3294nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h', 
3295                      'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
3296                      'R 3 c r','R -3 c r','R -3 m r',),
3297
3298#Use the space groups types in this order to list the symbols in the
3299#order they are listed in the International Tables, vol. A'''
3300symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 
3301               'trigonal', 'hexagonal', 'cubic')
3302
3303# self-test materials follow. Requires files in directory testinp
3304selftestlist = []
3305'''Defines a list of self-tests'''
3306selftestquiet = True
3307def _ReportTest():
3308    'Report name and doc string of current routine when ``selftestquiet`` is False'
3309    if not selftestquiet:
3310        import inspect
3311        caller = inspect.stack()[1][3]
3312        doc = eval(caller).__doc__
3313        if doc is not None:
3314            print('testing '+__file__+' with '+caller+' ('+doc+')')
3315        else:
3316            print('testing '+__file__()+" with "+caller)
3317def test0():
3318    '''self-test #0: exercise MoveToUnitCell'''
3319    _ReportTest()
3320    msg = "MoveToUnitCell failed"
3321    assert (MoveToUnitCell([1,2,3])[0] == [0,0,0]).all, msg
3322    assert (MoveToUnitCell([2,-1,-2])[0] == [0,0,0]).all, msg
3323    assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9)[0] < 1e-6, msg
3324    assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1)[0] < 1e-6, msg
3325selftestlist.append(test0)
3326
3327def test1():
3328    '''self-test #1: SpcGroup against previous results'''
3329    #'''self-test #1: SpcGroup and SGPrint against previous results'''
3330    _ReportTest()
3331    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3332    if ospath.exists(testdir):
3333        if testdir not in sys.path: sys.path.insert(0,testdir)
3334    import spctestinp
3335    def CompareSpcGroup(spc, referr, refdict, reflist): 
3336        'Compare output from GSASIIspc.SpcGroup with results from a previous run'
3337        # if an error is reported, the dictionary can be ignored
3338        msg0 = "CompareSpcGroup failed on space group %s" % spc
3339        result = SpcGroup(spc)
3340        if result[0] == referr and referr > 0: return True
3341#        #print result[1]['SpGrp']
3342        #msg = msg0 + " in list lengths"
3343        #assert len(keys) == len(refdict.keys()), msg
3344        for key in refdict.keys():
3345            if key == 'SGOps' or  key == 'SGCen':
3346                msg = msg0 + (" in key %s length" % key)
3347                assert len(refdict[key]) == len(result[1][key]), msg
3348                for i in range(len(refdict[key])):
3349                    msg = msg0 + (" in key %s level 0" % key)
3350                    assert np.allclose(result[1][key][i][0],refdict[key][i][0]), msg
3351                    msg = msg0 + (" in key %s level 1" % key)
3352                    assert np.allclose(result[1][key][i][1],refdict[key][i][1]), msg
3353            else:
3354                msg = msg0 + (" in key %s" % key)
3355                assert result[1][key] == refdict[key], msg
3356        msg = msg0 + (" in key %s reflist" % key)
3357        #for (l1,l2) in zip(reflist, SGPrint(result[1])):
3358        #    assert l2.replace('\t','').replace(' ','') == l1.replace(' ',''), 'SGPrint ' +msg
3359        # for now disable SGPrint testing, output has changed
3360        #assert reflist == SGPrint(result[1]), 'SGPrint ' +msg
3361    for spc in spctestinp.SGdat:
3362        CompareSpcGroup(spc, 0, spctestinp.SGdat[spc], spctestinp.SGlist[spc] )
3363selftestlist.append(test1)
3364
3365def test2():
3366    '''self-test #2: SpcGroup against cctbx (sgtbx) computations'''
3367    _ReportTest()
3368    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3369    if ospath.exists(testdir):
3370        if testdir not in sys.path: sys.path.insert(0,testdir)
3371    import sgtbxtestinp
3372    def CompareWcctbx(spcname, cctbx_in, debug=0):
3373        'Compare output from GSASIIspc.SpcGroup with results from cctbx.sgtbx'
3374        cctbx = cctbx_in[:] # make copy so we don't delete from the original
3375        spc = (SpcGroup(spcname))[1]
3376        if debug: print (spc['SpGrp'])
3377        if debug: print (spc['SGCen'])
3378        latticetype = spcname.strip().upper()[0]
3379        # lattice type of R implies Hexagonal centering", fix the rhombohedral settings
3380        if latticetype == "R" and len(spc['SGCen']) == 1: latticetype = 'P'
3381        assert latticetype == spc['SGLatt'], "Failed: %s does not match Lattice: %s" % (spcname, spc['SGLatt'])
3382        onebar = [1]
3383        if spc['SGInv']: onebar.append(-1)
3384        for (op,off) in spc['SGOps']:
3385            for inv in onebar:
3386                for cen in spc['SGCen']:
3387                    noff = off + cen
3388                    noff = MoveToUnitCell(noff)[0]
3389                    mult = tuple((op*inv).ravel().tolist())
3390                    if debug: print ("\n%s: %s + %s" % (spcname,mult,noff))
3391                    for refop in cctbx:
3392                        if debug: print (refop)
3393                        # check the transform
3394                        if refop[:9] != mult: continue
3395                        if debug: print ("mult match")
3396                        # check the translation
3397                        reftrans = list(refop[-3:])
3398                        reftrans = MoveToUnitCell(reftrans)[0]
3399                        if all(abs(noff - reftrans) < 1.e-5):
3400                            cctbx.remove(refop)
3401                            break
3402                    else:
3403                        assert False, "failed on %s:\n\t %s + %s" % (spcname,mult,noff)
3404    for key in sgtbxtestinp.sgtbx:
3405        CompareWcctbx(key, sgtbxtestinp.sgtbx[key])
3406selftestlist.append(test2)
3407
3408def test3(): 
3409    '''self-test #3: exercise SytSym (includes GetOprPtrName, GenAtom, GetKNsym)
3410     for selected space groups against info in IT Volume A '''
3411    _ReportTest()
3412    def ExerciseSiteSym (spc, crdlist):
3413        'compare site symmetries and multiplicities for a specified space group'
3414        msg = "failed on site sym test for %s" % spc
3415        (E,S) = SpcGroup(spc)
3416        assert not E, msg
3417        for t in crdlist:
3418            symb, m, n, od = SytSym(t[0],S)
3419            if symb.strip() != t[2].strip() or m != t[1]:
3420                print (spc,t[0],m,n,symb,t[2],od)
3421            assert m == t[1]
3422            #assert symb.strip() == t[2].strip()
3423
3424    ExerciseSiteSym('p 1',[
3425            ((0.13,0.22,0.31),1,'1'),
3426            ((0,0,0),1,'1'),
3427            ])
3428    ExerciseSiteSym('p -1',[
3429            ((0.13,0.22,0.31),2,'1'),
3430            ((0,0.5,0),1,'-1'),
3431            ])
3432    ExerciseSiteSym('C 2/c',[
3433            ((0.13,0.22,0.31),8,'1'),
3434            ((0.0,.31,0.25),4,'2(y)'),
3435            ((0.25,.25,0.5),4,'-1'),
3436            ((0,0.5,0),4,'-1'),
3437            ])
3438    ExerciseSiteSym('p 2 2 2',[
3439            ((0.13,0.22,0.31),4,'1'),
3440            ((0,0.5,.31),2,'2(z)'),
3441            ((0.5,.31,0.5),2,'2(y)'),
3442            ((.11,0,0),2,'2(x)'),
3443            ((0,0.5,0),1,'222'),
3444            ])
3445    ExerciseSiteSym('p 4/n',[
3446            ((0.13,0.22,0.31),8,'1'),
3447            ((0.25,0.75,.31),4,'2(z)'),
3448            ((0.5,0.5,0.5),4,'-1'),
3449            ((0,0.5,0),4,'-1'),
3450            ((0.25,0.25,.31),2,'4(001)'),
3451            ((0.25,.75,0.5),2,'-4(001)'),
3452            ((0.25,.75,0.0),2,'-4(001)'),
3453            ])
3454    ExerciseSiteSym('p 31 2 1',[
3455            ((0.13,0.22,0.31),6,'1'),
3456            ((0.13,0.0,0.833333333),3,'2(100)'),
3457            ((0.13,0.13,0.),3,'2(110)'),
3458            ])
3459    ExerciseSiteSym('R 3 c',[
3460            ((0.13,0.22,0.31),18,'1'),
3461            ((0.0,0.0,0.31),6,'3'),
3462            ])
3463    ExerciseSiteSym('R 3 c R',[
3464            ((0.13,0.22,0.31),6,'1'),
3465            ((0.31,0.31,0.31),2,'3(111)'),
3466            ])
3467    ExerciseSiteSym('P 63 m c',[
3468            ((0.13,0.22,0.31),12,'1'),
3469            ((0.11,0.22,0.31),6,'m(100)'),
3470            ((0.333333,0.6666667,0.31),2,'3m(100)'),
3471            ((0,0,0.31),2,'3m(100)'),
3472            ])
3473    ExerciseSiteSym('I a -3',[
3474            ((0.13,0.22,0.31),48,'1'),
3475            ((0.11,0,0.25),24,'2(x)'),
3476            ((0.11,0.11,0.11),16,'3(111)'),
3477            ((0,0,0),8,'-3(111)'),
3478            ])
3479selftestlist.append(test3)
3480
3481if __name__ == '__main__':
3482    # run self-tests
3483    selftestquiet = False
3484    for test in selftestlist:
3485        test()
3486    print ("OK")
Note: See TracBrowser for help on using the repository browser.