source: trunk/GSASIIspc.py @ 3765

Last change on this file since 3765 was 3765, checked in by vondreele, 3 years ago

allow a selected data copy of 'newLeBail' flag - facilitates sequential LeBail?
remove SawTooth? as modulation fxn - duplicates ZigZag?

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