source: trunk/GSASIIspc.py @ 3172

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

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