source: trunk/GSASIIspc.py @ 3171

Last change on this file since 3171 was 3171, checked in by vondreele, 6 years ago

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