source: trunk/GSASIIspc.py @ 3768

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

name fix for single crystal 1D plots
remove ZigZag? & Block from modulation derivatives - now done numerically
fix ZigZag? modulation function & proper limits for Tmin & Tmax in GUI

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 173.6 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-19 18:50:45 +0000 (Wed, 19 Dec 2018) $
12# $Author: vondreele $
13# $Revision: 3768 $
14# $URL: trunk/GSASIIspc.py $
15# $Id: GSASIIspc.py 3768 2018-12-19 18:50:45Z 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: 3768 $")
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(0. < t-(Tmm[0])%1. <= DT,-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        delt5 = np.ones(5)*0.001
2657        delt6 = np.eye(6)*0.001
2658        if 'Fourier' in waveType:
2659            dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2660              #3x6x12 modulated position array (X,Spos,tau)& force positive
2661        elif waveType in ['ZigZag','Block']:
2662            if waveType == 'ZigZag':
2663                dX = posZigZag(tau,delt5[:2],delt5[2:])
2664            else:
2665                dX = posBlock(tau,delt5[:2],delt5[2:])
2666        dXTP = []
2667        if siteSym == '1':
2668            CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2669        elif siteSym == '-1':
2670            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.],]
2671        else:
2672            if 'Fourier' in waveType:
2673                CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2674            elif waveType in ['ZigZag','Block']:
2675                CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2676                    np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2677            XSC = np.ones(6,dtype='i')
2678            dXTP = []
2679            for i in SdIndx:
2680                sop = Sop[i]
2681                ssop = SSop[i]
2682                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2683                xsc = np.ones(6,dtype='i')
2684                if 'Fourier' in waveType:
2685                    dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,nxs,nxs]
2686                elif waveType == 'ZigZag':
2687                    dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2688                elif waveType == 'Block':
2689                    dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2690                dXT = np.inner(sop[0],dXT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2691                dXT = np.swapaxes(dXT,1,2)      # back to array(3x6x49)
2692                dXT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2693                dXTP.append(dXT)
2694                if waveType == 'Fourier':
2695                    for i in range(3):
2696                        if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2697                            xsc[i] = 0
2698                        if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2699                            xsc[i+3] = 0
2700                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2701                        xsc[3:6] = 0
2702                        CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2703                            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2704                        if dT:
2705                            if '(x)' in siteSym:
2706                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2707                                if 'm' in siteSym and len(SdIndx) == 1:
2708                                    CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2709                            elif '(y)' in siteSym:
2710                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2711                                if 'm' in siteSym and len(SdIndx) == 1:
2712                                    CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2713                            elif '(z)' in siteSym:
2714                                CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2715                                if 'm' in siteSym and len(SdIndx) == 1:
2716                                    CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2717                        else:
2718                            CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2719                    if '4/mmm' in laue:
2720                        if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2721                            if '(xy)' in siteSym:
2722                                CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2723                                if dT:
2724                                    CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2725                                else:
2726                                    CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2727                        if '(xy)' in siteSym or '(+-0)' in siteSym:
2728                            mul = 1
2729                            if '(+-0)' in siteSym:
2730                                mul = -1
2731                            if np.allclose(dX[0,0,:],dXT[1,0,:]):
2732                                CSI[0][3:5] = [[11,0,0],[11,0,0]]
2733                                CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2734                                xsc[3:5] = 0
2735                            if np.allclose(dX[0,3,:],dXT[0,4,:]):
2736                                CSI[0][:2] = [[12,0,0],[12,0,0]]
2737                                CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2738                                xsc[:2] = 0
2739                else:
2740                    for i in range(3):
2741                        if not np.allclose(dX[:,i],dXT[i,:,i]):
2742                            xsc[i] = 0
2743                XSC &= xsc
2744                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2745            if waveType == 'Fourier':
2746                n = -1
2747                if debug: print (XSC)
2748                for i,X in enumerate(XSC):
2749                    if X:
2750                        n += 1
2751                        CSI[0][i][0] = n+1
2752                        CSI[1][i][0] = 1.0
2753           
2754        return list(CSI),dX,dXTP
2755       
2756    def DoUij():
2757        delt12 = np.eye(12)*0.0001
2758        dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
2759        dUTP = []
2760        if siteSym == '1':
2761            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2762                [7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2763        elif siteSym == '-1':
2764            CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],   \
2765                6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2766        else:
2767            CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2768            USC = np.ones(12,dtype='i')
2769            dUTP = []
2770            dtau = 0.
2771            for i in SdIndx:
2772                sop = Sop[i]
2773                ssop = SSop[i]
2774                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2775                usc = np.ones(12,dtype='i')
2776                dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
2777                dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3)    #convert dUT to 12x49x3x3
2778                dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2779                dUT = np.array(U2Uij(dUijT))    #convert to 6x12x49
2780                dUT = dUT[:,:,np.argsort(tauT)]
2781                dUT[:,:6,:] *=(ssdet*sdet)
2782                dUTP.append(dUT)
2783                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2784                    if dT:
2785                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2786                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2787                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2788                        [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2789                    else:
2790                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2791                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2792                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2793                        [0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2794                    if 'mm2(x)' in siteSym and dT:
2795                        CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2796                        USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2797                    elif '(xy)' in siteSym and dT:
2798                        CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2799                            [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2800                        CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2801                        USC = [1,1,1,1,1,1,1,1,1,1,1,1]                             
2802                    elif '(x)' in siteSym and dT:
2803                        CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2804                    elif '(y)' in siteSym and dT:
2805                        CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2806                    elif '(z)' in siteSym and dT:
2807                        CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2808                    for i in range(6):
2809                        if not USC[i]:
2810                            CSI[0][i] = [0,0,0]
2811                            CSI[1][i] = [0.,0.,0.]
2812                            CSI[0][i+6] = [0,0,0]
2813                            CSI[1][i+6] = [0.,0.,0.]
2814                else:                       
2815                    for i in range(6):
2816                        if not np.allclose(dU[i,i,:],dUT[i,i,:]):  #sin part
2817                            usc[i] = 0
2818                        if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]):   #cos part
2819                            usc[i+6] = 0
2820                    if np.any(dUT[1,0,:]):
2821                        if '4/m' in siteSym:
2822                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2823                            if ssop[1][3]:
2824                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2825                                usc[9] = 1
2826                            else:
2827                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2828                                usc[9] = 0
2829                        elif '4' in siteSym:
2830                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2831                            CSI[0][:2] = [[11,0,0],[11,0,0]]
2832                            if ssop[1][3]:
2833                                CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2834                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2835                                usc[2] = 0
2836                                usc[8] = 0
2837                                usc[3] = 1
2838                                usc[9] = 1
2839                            else:
2840                                CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2841                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2842                                usc[2] = 1
2843                                usc[8] = 1
2844                                usc[3] = 0               
2845                                usc[9] = 0
2846                        elif 'xy' in siteSym or '+-0' in siteSym:
2847                            if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2848                                CSI[0][4:6] = [[12,0,0],[12,0,0]]
2849                                CSI[0][6:8] = [[11,0,0],[11,0,0]]
2850                                CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2851                                CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2852                                usc[4:6] = 0
2853                                usc[6:8] = 0
2854                           
2855                    if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2856                USC &= usc
2857            if debug: print (USC)
2858            if not np.any(dtau%.5):
2859                n = -1
2860                for i,U in enumerate(USC):
2861                    if U:
2862                        n += 1
2863                        CSI[0][i][0] = n+1
2864                        CSI[1][i][0] = 1.0
2865   
2866        return list(CSI),dU,dUTP
2867   
2868    def DoMag():
2869        delt6 = np.eye(6)*0.001
2870        dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,nxs,nxs]
2871        dMTP = []
2872        CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2873        if siteSym == '1':
2874            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2875        elif siteSym in ['-1','mmm',]:
2876            CSI = 3*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3*[[0.,0.,0.],]+3*[[1.,0.,0.],]
2877        elif siteSym in ['4(z)','422(z)']:
2878            CSI[0][0][0] = CSI[0][4][1] = 1
2879            CSI[1][0][0] = 1.0
2880            CSI[1][4][1] = -1.0
2881        elif siteSym in ['-4m2(z)','422(z)',]:
2882            CSI[0][5][0] = 1
2883            CSI[1][5][0] = 1.0
2884        elif siteSym in ['-32(100)','-3',]:
2885            CSI[0][2][0] = 1
2886            CSI[1][2][0] = 1.0
2887        elif siteSym in ['3',]:
2888            CSI[0][0][0] = CSI[0][3][0] = CSI[0][4][0] = 1
2889            CSI[1][0][0] = -np.sqrt(3.0)
2890            CSI[1][3][0] = 2.0
2891            CSI[1][4][0] = 1.0
2892        elif siteSym in ['622','2(100)','32(100)',]:
2893            CSI[0][0][0] = CSI[0][1][0] = CSI[0][3][0] = 1
2894            CSI[1][0][0] = 1.0
2895            CSI[1][1][0] = 2.0
2896            CSI[1][3][0] = np.sqrt(3.0)
2897        else:
2898              #3x6x12 modulated moment array (M,Spos,tau)& force positive
2899            CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2900            MSC = np.ones(6,dtype='i')
2901            dMTP = []
2902            for i in SdIndx:
2903                sop = Sop[i]
2904                ssop = SSop[i]
2905                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2906                msc = np.ones(6,dtype='i')
2907                dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,nxs,nxs]
2908                dMT = np.inner(sop[0],dMT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2909                dMT = np.swapaxes(dMT,1,2)      # back to array(3x6x49)
2910                dMT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2911                dMTP.append(dMT)
2912                for i in range(3):
2913                    if not np.allclose(dM[i,i,:],sdet*dMT[i,i,:]):
2914                        msc[i] = 0
2915                    if not np.allclose(dM[i,i+3,:],sdet*dMT[i,i+3,:]):
2916                        msc[i+3] = 0
2917                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2918                    msc[3:6] = 0
2919                    CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2920                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2921                    if dT:
2922                        if '(x)' in siteSym:
2923                            CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2924                            if 'm' in siteSym and len(SdIndx) == 1:
2925                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2926                        elif '(y)' in siteSym:
2927                            CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2928                            if 'm' in siteSym and len(SdIndx) == 1:
2929                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2930                        elif '(z)' in siteSym:
2931                            CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2932                            if 'm' in siteSym and len(SdIndx) == 1:
2933                                CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2934                    else:
2935                        CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2936                if '4/mmm' in laue:
2937                    if siteSym in ['4/mmm(z)',]:
2938                        CSI = 3*[[0,0,0],]+[[0,0,0],[0,0,0],[1,0,0]],3*[[0.,0.,0.],]+3*[[1.,0.,0.],]
2939                    if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2940                        if '(xy)' in siteSym:
2941                            CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2942                            if dT:
2943                                CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2944                            else:
2945                                CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2946                    if '(xy)' in siteSym or '(+-0)' in siteSym:
2947                        mul = 1
2948                        if '(+-0)' in siteSym:
2949                            mul = -1
2950                        if np.allclose(dM[0,0,:],dMT[1,0,:]):
2951                            CSI[0][3:5] = [[11,0,0],[11,0,0]]
2952                            CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2953                            msc[3:5] = 0
2954                        if np.allclose(dM[0,3,:],dMT[0,4,:]):
2955                            CSI[0][:2] = [[12,0,0],[12,0,0]]
2956                            CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2957                            msc[:2] = 0
2958                MSC &= msc
2959                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2960            n = -1
2961            if debug: print (MSC)
2962            for i,M in enumerate(MSC):
2963                if M:
2964                    n += 1
2965                    CSI[0][i][0] = n+1
2966                    CSI[1][i][0] = 1.0
2967
2968        return list(CSI),dM,dMTP
2969       
2970    if debug: print ('super space group: '+SSGData['SSpGrp'])
2971    xyz = np.array(XYZ)%1.
2972    SGOps = copy.deepcopy(SGData['SGOps'])
2973    laue = SGData['SGLaue']
2974    siteSym = SytSym(XYZ,SGData)[0].strip()
2975    if debug: print ('siteSym: '+siteSym)
2976    SSGOps = copy.deepcopy(SSGData['SSGOps'])
2977    #expand ops to include inversions if any
2978    if SGData['SGInv'] and not SGData['SGFixed']:
2979        for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
2980            SGOps.append([-op[0],-op[1]%1.])
2981            SSGOps.append([-sop[0],-sop[1]%1.])
2982    #build set of sym ops around special position       
2983    SSop = []
2984    Sop = []
2985    Sdtau = []
2986    for iop,Op in enumerate(SGOps):         
2987        nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
2988        if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
2989            SSop.append(SSGOps[iop])
2990            Sop.append(SGOps[iop])
2991            ssopinv = nl.inv(SSGOps[iop][0])
2992            mst = ssopinv[3][:3]
2993            epsinv = ssopinv[3][3]
2994            Sdtau.append(np.sum(mst*(XYZ-SGOps[iop][1])-epsinv*SSGOps[iop][1][3]))
2995    SdIndx = np.argsort(np.array(Sdtau))     # just to do in sensible order
2996    if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
2997    #setup displacement arrays
2998    tau = np.linspace(-1,1,49,True)
2999    #make modulation arrays - one parameter at a time
3000    if Stype == 'Sfrac':
3001        CSI,dF,dFTP = DoFrac()
3002    elif Stype == 'Spos':
3003        CSI,dF,dFTP = DoXYZ()
3004        CSI[0] = orderParms(CSI[0])
3005    elif Stype == 'Sadp':
3006        CSI,dF,dFTP = DoUij()
3007        CSI[0] = orderParms(CSI[0]) 
3008    elif Stype == 'Smag':
3009        CSI,dF,dFTP = DoMag()
3010
3011    if debug:
3012        return CSI,dF,dFTP
3013    else:
3014        return CSI,[],[]
3015   
3016def MustrainNames(SGData):
3017    'Needs a doc string'
3018    laue = SGData['SGLaue']
3019    uniq = SGData['SGUniq']
3020    if laue in ['m3','m3m']:
3021        return ['S400','S220']
3022    elif laue in ['6/m','6/mmm','3m1']:
3023        return ['S400','S004','S202']
3024    elif laue in ['31m','3']:
3025        return ['S400','S004','S202','S211']
3026    elif laue in ['3R','3mR']:
3027        return ['S400','S220','S310','S211']
3028    elif laue in ['4/m','4/mmm']:
3029        return ['S400','S004','S220','S022']
3030    elif laue in ['mmm']:
3031        return ['S400','S040','S004','S220','S202','S022']
3032    elif laue in ['2/m']:
3033        SHKL = ['S400','S040','S004','S220','S202','S022']
3034        if uniq == 'a':
3035            SHKL += ['S013','S031','S211']
3036        elif uniq == 'b':
3037            SHKL += ['S301','S103','S121']
3038        elif uniq == 'c':
3039            SHKL += ['S130','S310','S112']
3040        return SHKL
3041    else:
3042        SHKL = ['S400','S040','S004','S220','S202','S022']
3043        SHKL += ['S310','S103','S031','S130','S301','S013']
3044        SHKL += ['S211','S121','S112']
3045        return SHKL
3046       
3047def HStrainVals(HSvals,SGData):
3048    laue = SGData['SGLaue']
3049    uniq = SGData['SGUniq']
3050    DIJ = np.zeros(6)
3051    if laue in ['m3','m3m']:
3052        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
3053    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3054        DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
3055    elif laue in ['3R','3mR']:
3056        DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
3057    elif laue in ['4/m','4/mmm']:
3058        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
3059    elif laue in ['mmm']:
3060        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3061    elif laue in ['2/m']:
3062        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3063        if uniq == 'a':
3064            DIJ[5] = HSvals[3]
3065        elif uniq == 'b':
3066            DIJ[4] = HSvals[3]
3067        elif uniq == 'c':
3068            DIJ[3] = HSvals[3]
3069    else:
3070        DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
3071    return DIJ
3072
3073def HStrainNames(SGData):
3074    'Needs a doc string'
3075    laue = SGData['SGLaue']
3076    uniq = SGData['SGUniq']
3077    if laue in ['m3','m3m']:
3078        return ['D11','eA']         #add cubic strain term
3079    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3080        return ['D11','D33']
3081    elif laue in ['3R','3mR']:
3082        return ['D11','D12']
3083    elif laue in ['4/m','4/mmm']:
3084        return ['D11','D33']
3085    elif laue in ['mmm']:
3086        return ['D11','D22','D33']
3087    elif laue in ['2/m']:
3088        Dij = ['D11','D22','D33']
3089        if uniq == 'a':
3090            Dij += ['D23']
3091        elif uniq == 'b':
3092            Dij += ['D13']
3093        elif uniq == 'c':
3094            Dij += ['D12']
3095        return Dij
3096    else:
3097        Dij = ['D11','D22','D33','D12','D13','D23']
3098        return Dij
3099   
3100def MustrainCoeff(HKL,SGData):
3101    'Needs a doc string'
3102    #NB: order of terms is the same as returned by MustrainNames
3103    laue = SGData['SGLaue']
3104    uniq = SGData['SGUniq']
3105    h,k,l = HKL
3106    Strm = []
3107    if laue in ['m3','m3m']:
3108        Strm.append(h**4+k**4+l**4)
3109        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
3110    elif laue in ['6/m','6/mmm','3m1']:
3111        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
3112        Strm.append(l**4)
3113        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
3114    elif laue in ['31m','3']:
3115        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
3116        Strm.append(l**4)
3117        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
3118        Strm.append(4.0*h*k*l*(h+k))
3119    elif laue in ['3R','3mR']:
3120        Strm.append(h**4+k**4+l**4)
3121        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
3122        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
3123        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
3124    elif laue in ['4/m','4/mmm']:
3125        Strm.append(h**4+k**4)
3126        Strm.append(l**4)
3127        Strm.append(3.0*(h*k)**2)
3128        Strm.append(3.0*((h*l)**2+(k*l)**2))
3129    elif laue in ['mmm']:
3130        Strm.append(h**4)
3131        Strm.append(k**4)
3132        Strm.append(l**4)
3133        Strm.append(3.0*(h*k)**2)
3134        Strm.append(3.0*(h*l)**2)
3135        Strm.append(3.0*(k*l)**2)
3136    elif laue in ['2/m']:
3137        Strm.append(h**4)
3138        Strm.append(k**4)
3139        Strm.append(l**4)
3140        Strm.append(3.0*(h*k)**2)
3141        Strm.append(3.0*(h*l)**2)
3142        Strm.append(3.0*(k*l)**2)
3143        if uniq == 'a':
3144            Strm.append(2.0*k*l**3)
3145            Strm.append(2.0*l*k**3)
3146            Strm.append(4.0*k*l*h**2)
3147        elif uniq == 'b':
3148            Strm.append(2.0*l*h**3)
3149            Strm.append(2.0*h*l**3)
3150            Strm.append(4.0*h*l*k**2)
3151        elif uniq == 'c':
3152            Strm.append(2.0*h*k**3)
3153            Strm.append(2.0*k*h**3)
3154            Strm.append(4.0*h*k*l**2)
3155    else:
3156        Strm.append(h**4)
3157        Strm.append(k**4)
3158        Strm.append(l**4)
3159        Strm.append(3.0*(h*k)**2)
3160        Strm.append(3.0*(h*l)**2)
3161        Strm.append(3.0*(k*l)**2)
3162        Strm.append(2.0*k*h**3)
3163        Strm.append(2.0*h*l**3)
3164        Strm.append(2.0*l*k**3)
3165        Strm.append(2.0*h*k**3)
3166        Strm.append(2.0*l*h**3)
3167        Strm.append(2.0*k*l**3)
3168        Strm.append(4.0*k*l*h**2)
3169        Strm.append(4.0*h*l*k**2)
3170        Strm.append(4.0*k*h*l**2)
3171    return Strm
3172
3173def MuShklMean(SGData,Amat,Shkl):
3174   
3175    def genMustrain(xyz,Shkl):
3176        uvw = np.inner(Amat.T,xyz)
3177        Strm = np.array(MustrainCoeff(uvw,SGData))
3178        Sum = np.sum(np.multiply(Shkl,Strm))
3179        Sum = np.where(Sum > 0.01,Sum,0.01)
3180        Sum = np.sqrt(Sum)
3181        return Sum*xyz
3182       
3183    PHI = np.linspace(0.,360.,30,True)
3184    PSI = np.linspace(0.,180.,30,True)
3185    X = np.outer(npcosd(PHI),npsind(PSI))
3186    Y = np.outer(npsind(PHI),npsind(PSI))
3187    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3188    XYZ = np.dstack((X,Y,Z))
3189    XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
3190    return np.sqrt(np.sum(XYZ**2)/900.)
3191   
3192def Muiso2Shkl(muiso,SGData,cell):
3193    "this is to convert isotropic mustrain to generalized Shkls"
3194    import GSASIIlattice as G2lat
3195    A = G2lat.cell2AB(cell)[0]
3196   
3197    def minMus(Shkl,muiso,H,SGData,A):
3198        U = np.inner(A.T,H)
3199        S = np.array(MustrainCoeff(U,SGData))
3200        Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,nxs]),axis=0))
3201        rad = np.sqrt(np.sum((Sum[:,nxs]*H)**2,axis=1))
3202        return (muiso-rad)**2
3203       
3204    laue = SGData['SGLaue']
3205    PHI = np.linspace(0.,360.,60,True)
3206    PSI = np.linspace(0.,180.,60,True)
3207    X = np.outer(npsind(PHI),npsind(PSI))
3208    Y = np.outer(npcosd(PHI),npsind(PSI))
3209    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3210    HKL = np.dstack((X,Y,Z))
3211    if laue in ['m3','m3m']:
3212        S0 = [1000.,1000.]
3213    elif laue in ['6/m','6/mmm','3m1']:
3214        S0 = [1000.,1000.,1000.]
3215    elif laue in ['31m','3']:
3216        S0 = [1000.,1000.,1000.,1000.]
3217    elif laue in ['3R','3mR']:
3218        S0 = [1000.,1000.,1000.,1000.]
3219    elif laue in ['4/m','4/mmm']:
3220        S0 = [1000.,1000.,1000.,1000.]
3221    elif laue in ['mmm']:
3222        S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
3223    elif laue in ['2/m']:
3224        S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
3225    else:
3226        S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000., 
3227            1000.,1000.,0.,0.,0.]
3228    S0 = np.array(S0)
3229    HKL = np.reshape(HKL,(-1,3))
3230    result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
3231    return result[0]
3232       
3233def PackRot(SGOps):
3234    IRT = []
3235    for ops in SGOps:
3236        M = ops[0]
3237        irt = 0
3238        for j in range(2,-1,-1):
3239            for k in range(2,-1,-1):
3240                irt *= 3
3241                irt += M[k][j]
3242        IRT.append(int(irt))
3243    return IRT
3244       
3245def SytSym(XYZ,SGData):
3246    '''
3247    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
3248
3249    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
3250    :param SGData: from SpcGroup
3251    :Returns: a four element tuple:
3252
3253     * The 1st element is a code for the site symmetry (see GetKNsym)
3254     * The 2nd element is the site multiplicity
3255     * Ndup number of overlapping operators
3256     * dupDir Dict - dictionary of overlapping operators
3257
3258    '''
3259    Mult = 1
3260    Isym = 0
3261    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
3262        Isym = 1073741824
3263    Jdup = 0
3264    Ndup = 0
3265    dupDir = {}
3266    inv = SGData['SGInv']+1
3267    icen = SGData['SGCen']
3268    Ncen = len(icen)
3269    if SGData['SGFixed']:       #already in list of operators
3270        inv = 1
3271    if SGData['SGGray'] and Ncen > 1: Ncen //= 2
3272    Xeqv = list(GenAtom(XYZ,SGData,True))
3273#    for xeqv in Xeqv:   print(xeqv)
3274    IRT = PackRot(SGData['SGOps'])
3275    L = -1
3276    for ic,cen in enumerate(icen[:Ncen]):
3277        for invers in range(int(inv)):
3278            for io,ops in enumerate(SGData['SGOps']):
3279                irtx = (1-2*invers)*IRT[io]
3280                L += 1
3281                if not Xeqv[L][1]:
3282                    Ndup = io
3283                    Jdup += 1
3284                    jx = GetOprPtrNumber(str(irtx))   #[KN table no,op name,KNsym ptr]
3285                    if jx < 39:
3286                        px = GetOprName(str(irtx))
3287                        if Xeqv[L][-1] < 0:
3288                            if '(' in px:
3289                                px = px.split('(')
3290                                px[0] += "'"
3291                                px = '('.join(px)
3292                            else:   
3293                                px += "'"
3294                        dupDir[px] = L
3295                        Isym += 2**(jx-1)
3296    if Isym == 1073741824: Isym = 0
3297    Mult = len(SGData['SGOps'])*Ncen*inv//Jdup
3298         
3299    return GetKNsym(str(Isym)),Mult,Ndup,dupDir
3300   
3301def MagSytSym(SytSym,dupDir,SGData):
3302    '''
3303    site sym operations: 1,-1,2,3,-3,4,-4,6,-6,m need to be marked if spin inversion
3304    '''
3305    SGData['GenSym'],SGData['GenFlg'] = GetGenSym(SGData)[:2]
3306#    print('SGPtGrp',SGData['SGPtGrp'],'SytSym',SytSym,'MagSpGrp',SGData['MagSpGrp'])
3307#    print('dupDir',dupDir)
3308    SplitSytSym = SytSym.split('(')
3309    if SGData['SGGray']:
3310        return SytSym+"1'"
3311    if SytSym == '1':       #genersl position
3312        return SytSym
3313    if SplitSytSym[0] == SGData['SGPtGrp']:     #simple cases
3314        try:
3315            MagSytSym = SGData['MagSpGrp'].split()[1]
3316        except IndexError:
3317            MagSytSym = SGData['MagSpGrp'][1:].strip("1'")
3318        if len(SplitSytSym) > 1:
3319            MagSytSym += '('+SplitSytSym[1]
3320        return MagSytSym
3321    if len(dupDir) == 1:
3322        return list(dupDir.keys())[0]
3323   
3324   
3325    if '2/m' in SytSym:         #done I think; last 2wo might be not needed
3326        ops = {'(x)':['2(x)','m(x)'],'(y)':['2(y)','m(y)'],'(z)':['2(z)','m(z)'],
3327               '(100)':['2(100)','m(100)'],'(010)':['2(010)','m(010)'],'(001)':['2(001)','m(001)'],
3328               '(120)':['2(120)','m(120)'],'(210)':['2(210)','m(210)'],'(+-0)':['2(+-0)','m(+-0)'],
3329               '(110)':['2(110)','m(110)']}
3330   
3331    elif '4/mmm' in SytSym:
3332        ops = {'(x)':['4(x)','m(x)','m(y)','m(0+-)'],   #m(0+-) for cubic m3m?
3333               '(y)':['4(y)','m(y)','m(z)','m(+0-)'],   #m(+0-)
3334               '(z)':['4(z)','m(z)','m(x)','m(+-0)']}   #m(+-0)
3335    elif '4mm' in SytSym:
3336        ops = {'(x)':['4(x)','m(y)','m(yz)'],'(y)':['4(y)','m(z)','m(xz)'],'(z)':['4(z)','m(x)','m(110)']}
3337    elif '422' in SytSym:
3338        ops = {'(x)':['4(x)','2(y)','2(yz)'],'(y)':['4(y)','2(z)','2(xz)'],'(z)':['4(z)','2(x)','2(110)']}
3339    elif '-4m2' in SytSym:
3340        ops = {'(x)':['-4(x)','m(x)','2(yz)'],'(y)':['-4(y)','m(y)','2(xz)'],'(z)':['-4(z)','m(z)','2(110)']}
3341    elif '-42m' in SytSym:
3342        ops = {'(x)':['-4(x)','2(y)','m(yz)'],'(y)':['-4(y)','2(z)','m(xz)'],'(z)':['-4(z)','2(x)','m(110)']}
3343    elif '-4' in SytSym:
3344        ops = {'(x)':['-4(x)',],'(y)':['-4(y)',],'(z)':['-4(z)',],}
3345    elif '4' in SytSym:
3346        ops = {'(x)':['4(x)',],'(y)':['4(y)',],'(z)':['4(z)',],}
3347
3348    elif '222' in SytSym:
3349        ops = {'':['2(x)','2(y)','2(z)'],
3350                   '(x)':['2(y)','2(z)','2(x)'],'(y)':['2(x)','2(z)','2(y)'],'(z)':['2(x)','2(y)','2(z)'],
3351                   '(100)':['2(z)','2(100)','2(120)',],'(010)':['2(z)','2(010)','2(210)',],
3352                   '(110)':['2(z)','2(110)','2(+-0)',],}
3353    elif 'mm2' in SytSym:
3354        ops = {'(x)':['m(y)','m(z)','2(x)'],'(y)':['m(x)','m(z)','2(y)'],'(z)':['m(x)','m(y)','2(z)'],
3355               '(xy)':['m(+-0)','m(z)','2(110)'],'(yz)':['m(0+-)','m(xz)','2(yz)'],     #not 2(xy)!
3356               '(xz)':['m(+0-)','m(y)','2(xz)'],'(z100)':['m(100)','m(120)','2(z)'],
3357               '(z010)':['m(010)','m(210)','2(z)'],'(z110)':['m(110)','m(+-0)','2(z)'],
3358               '(+-0)':[ 'm(110)','m(z)','2(+-0)'],'(d100)':['m(yz)','m(0+-)','2(xz)'],
3359               '(d010)':['m(xz)','m(+0-)','2(y)'],'(d001)':['m(110)','m(+-0)','2(z)'],
3360               '(210)':['m(z)','m(010)','2(210)'],'(120)':['m(z)','m(100)','2(120)'],
3361               '(100)':['m(z)','m(120)','2(100)',],'(010)':['m(z)','m(210)','2(010)',],
3362               '(110)':['m(z)','m(+-0)','2(110)',],}
3363    elif 'mmm' in SytSym:
3364        ops = {'':['m(x)','m(y)','m(z)'],
3365                   '(100)':['m(z)','m(100)','m(120)',],'(010)':['m(z)','m(010)','m(210)',],
3366                   '(110)':['m(z)','m(110)','m(+-0)',],
3367                   '(x)':['m(x)','m(y)','m(z)'],'(y)':['m(x)','m(y)','m(z)'],'(z)':['m(x)','m(y)','m(z)'],}
3368       
3369    elif '32' in SytSym:
3370        ops = {'(120)':['3','2(120)',],'(100)':['3','2(100)'],'(111)':['3(111)','2(x)']}
3371    elif '23' in SytSym:
3372        ops = {'':['2(x)','3(111)']}
3373    elif 'm3' in SytSym:
3374        ops = {'(100)':['(+-0)',],'(+--)':[],'(-+-)':[],'(--+)':[]}
3375    elif '3m' in SytSym:
3376        ops = {'(111)':['3(111)','m(+-0)',],'(+--)':['3(+--)','m(0+-)',],
3377               '(-+-)':['3(-+-)','m(+0-)',],'(--+)':['3(--+)','m(+-0)',],
3378               '(100)':['3','m(100)'],'(120)':['3','m(210)',]}
3379   
3380    if SytSym.split('(')[0] in ['6/m','6mm','-6m2','622','-6','-3','-3m','-43m',]:     #not simple cases
3381        MagSytSym = SytSym
3382        if "-1'" in dupDir:
3383            if '-6' in SytSym:
3384                MagSytSym = MagSytSym.replace('-6',"-6'")
3385            elif '-3m' in SytSym:
3386                MagSytSym = MagSytSym.replace('-3m',"-3'm'")
3387            elif '-3' in SytSym:
3388                MagSytSym = MagSytSym.replace('-3',"-3'")
3389        elif '-6m2' in SytSym:
3390            if "m'(110)" in dupDir:
3391                MagSytSym = "-6m'2'("+SytSym.split('(')[1]
3392        elif '6/m' in SytSym:
3393            if "m'(z)" in dupDir:
3394                MagSytSym = "6'/m'"
3395        elif '6mm' in SytSym:
3396            if "m'(110)" in dupDir:
3397                MagSytSym = "6'm'm"
3398        elif '-43m' in SytSym:
3399            if "m'(110)" in dupDir:
3400                MagSytSym = "-43m'"
3401        return MagSytSym
3402    try:
3403        axis = '('+SytSym.split('(')[1]
3404    except IndexError:
3405        axis = ''
3406    MagSytSym = ''
3407    for m in ops[axis]:
3408        if m in dupDir:
3409            MagSytSym += m.split('(')[0]
3410        else:
3411            MagSytSym += m.split('(')[0]+"'"
3412        if '2/m' in SytSym and '2' in m:
3413            MagSytSym += '/'
3414        if '-3/m' in SytSym:
3415            MagSytSym = '-'+MagSytSym
3416       
3417    MagSytSym += axis
3418# some exceptions & special rules         
3419    if MagSytSym == "4'/m'm'm'": MagSytSym = "4/m'm'm'"
3420    return MagSytSym
3421   
3422#    if len(GenSym) == 3:
3423#        if SGSpin[1] < 0:
3424#            if 'mm2' in SytSym:
3425#                MagSytSym = "m'm'2"+'('+SplitSytSym[1]
3426#            else:   #bad rule for I41/a
3427#                MagSytSym = SplitSytSym[0]+"'"
3428#                if len(SplitSytSym) > 1:
3429#                    MagSytSym += '('+SplitSytSym[1]
3430#        else:
3431#            MagSytSym = SytSym
3432#        if len(SplitSytSym) >1:
3433#            if "-4'"+'('+SplitSytSym[1] in dupDir:
3434#                MagSytSym = MagSytSym.replace('-4',"-4'")
3435#            if "-6'"+'('+SplitSytSym[1] in dupDir:
3436#                MagSytSym = MagSytSym.replace('-6',"-6'")
3437#        return MagSytSym
3438#           
3439    return SytSym
3440
3441def UpdateSytSym(Phase):
3442    ''' Update site symmetry/site multiplicity after space group/BNS lattice change
3443    '''
3444    generalData = Phase['General']
3445    SGData = generalData['SGData']
3446    Atoms = Phase['Atoms']
3447    cx,ct,cs,cia = generalData['AtomPtrs']
3448    for atom in Atoms:
3449        XYZ = atom[cx:cx+3]
3450        sytsym,Mult = SytSym(XYZ,SGData)[:2]
3451        sytSym,Mul,Nop,dupDir = SytSym(XYZ,SGData)
3452        atom[cs] = sytsym
3453        if generalData['Type'] == 'magnetic':
3454            magSytSym = MagSytSym(sytSym,dupDir,SGData)
3455            atom[cs] = magSytSym
3456        atom[cs+1] = Mult
3457    return
3458   
3459def ElemPosition(SGData):
3460    ''' Under development.
3461    Object here is to return a list of symmetry element types and locations suitable
3462    for say drawing them.
3463    So far I have the element type... getting all possible locations without lookup may be impossible!
3464    '''
3465    Inv = SGData['SGInv']
3466    eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
3467    # get operators & expand if centrosymmetric
3468    SymElements = []
3469    Ops = SGData['SGOps']
3470    opM = np.array([op[0].T for op in Ops])
3471    opT = np.array([op[1] for op in Ops])
3472    if Inv:
3473        opM = np.concatenate((opM,-opM))
3474        opT = np.concatenate((opT,-opT))
3475    opMT = list(zip(opM,opT))
3476    for M,T in opMT[1:]:        #skip I
3477        Dt = int(nl.det(M))
3478        Tr = int(np.trace(M))
3479        Dt = -(Dt-1)//2
3480        Es = eleSym[Tr][Dt]
3481        if Dt:              #rotation-inversion
3482            I = np.eye(3)
3483            if Tr == 1:     #mirrors/glides
3484                if np.any(T):       #glide
3485                    M2 = np.inner(M,M)
3486                    MT = np.inner(M,T)+T
3487                    print ('glide',Es,MT)
3488                    print (M2)
3489                else:               #mirror
3490                    print ('mirror',Es,T)
3491                    print (I-M)
3492                X = [-1,-1,-1]
3493            elif Tr == -3:  # pure inversion
3494                X = np.inner(nl.inv(I-M),T)
3495                print ('inversion',Es,X)
3496            else:           #other rotation-inversion
3497                M2 = np.inner(M,M)
3498                MT = np.inner(M,T)+T
3499                print ('rot-inv',Es,MT)
3500                print (M2)
3501                X = [-1,-1,-1]
3502        else:               #rotations
3503            print ('rotation',Es)
3504            X = [-1,-1,-1]
3505        SymElements.append([Es,X])
3506       
3507    return SymElements
3508   
3509def ApplyStringOps(A,SGData,X,Uij=[]):
3510    'Needs a doc string'
3511    SGOps = SGData['SGOps']
3512    SGCen = SGData['SGCen']
3513    Ax = A.split('+')
3514    Ax[0] = int(Ax[0])
3515    iC = 1
3516    if Ax[0] < 0:
3517        iC = -1
3518    Ax[0] = abs(Ax[0])
3519    nA = Ax[0]%100-1
3520    cA = Ax[0]//100
3521    Cen = SGCen[cA]
3522    M,T = SGOps[nA]
3523    if len(Ax)>1:
3524        cellA = Ax[1].split(',')
3525        cellA = np.array([int(a) for a in cellA])
3526    else:
3527        cellA = np.zeros(3)
3528    newX = Cen+iC*(np.inner(M,X).T+T)+cellA
3529    if len(Uij):
3530        U = Uij2U(Uij)
3531        U = np.inner(M,np.inner(U,M).T)
3532        newUij = U2Uij(U)
3533        return [newX,newUij]
3534    else:
3535        return newX
3536       
3537def ApplyStringOpsMom(A,SGData,Mom):
3538    'Needs a doc string'
3539    SGOps = SGData['SGOps']
3540    Ax = A.split('+')
3541    Ax[0] = int(Ax[0])
3542    iAx = abs(Ax[0])
3543    nA = iAx%100-1
3544    nC = len(SGOps)*(iAx//100)
3545    NA = nA
3546    if Ax[0] < 0:
3547        NA += len(SGOps)
3548    M,T = SGOps[nA]
3549    if SGData['SGGray']:
3550        newMom = -np.inner(Mom,M).T*nl.det(M)*SGData['SpnFlp'][NA+nC]
3551    else:
3552        newMom = np.inner(Mom,M).T*nl.det(M)*SGData['SpnFlp'][NA+nC]
3553#        print(len(SGOps),Ax[0],iAx,nC,nA,NA,MT2text([M,T]).replace(' ',''),SGData['SpnFlp'][NA],Mom,newMom)
3554#    print(Mom,newMom,MT2text([M,T]),)
3555    return newMom
3556       
3557def StringOpsProd(A,B,SGData):
3558    """
3559    Find A*B where A & B are in strings '-' + '100*c+n' + '+ijk'
3560    where '-' indicates inversion, c(>0) is the cell centering operator,
3561    n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3562    Should return resultant string - C. SGData - dictionary using entries:
3563
3564       *  'SGCen': cell centering vectors [0,0,0] at least
3565       *  'SGOps': symmetry operations as [M,T] so that M*x+T = x'
3566
3567    """
3568    SGOps = SGData['SGOps']
3569    SGCen = SGData['SGCen']
3570    #1st split out the cell translation part & work on the operator parts
3571    Ax = A.split('+'); Bx = B.split('+')
3572    Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
3573    iC = 0
3574    if Ax[0]*Bx[0] < 0:
3575        iC = 1
3576    Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
3577    nA = Ax[0]%100-1;  nB = Bx[0]%100-1
3578    cA = Ax[0]//100;  cB = Bx[0]//100
3579    Cen = (SGCen[cA]+SGCen[cB])%1.0
3580    cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
3581    Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
3582    Mc = np.inner(Ma,Mb.T)
3583#    print Ma,Mb,Mc
3584    Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
3585#    print Ta,Tb,Tc
3586#    print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
3587    nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
3588    #now the cell translation part
3589    if len(Ax)>1:
3590        cellA = Ax[1].split(',')
3591        cellA = [int(a) for a in cellA]
3592    else:
3593        cellA = [0,0,0]
3594    if len(Bx)>1:
3595        cellB = Bx[1].split(',')
3596        cellB = [int(b) for b in cellB]
3597    else:
3598        cellB = [0,0,0]
3599    cellC = np.add(cellA,cellB)
3600    C = str(((nC+1)+(100*cC))*(1-2*iC))+'+'+ \
3601        str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
3602    return C
3603           
3604def U2Uij(U):
3605    #returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
3606    return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
3607   
3608def Uij2U(Uij):
3609    #returns the thermal motion tensor U from Uij as numpy array
3610    return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
3611
3612def StandardizeSpcName(spcgroup):
3613    '''Accept a spacegroup name where spaces may have not been used
3614    in the names according to the GSAS convention (spaces between symmetry
3615    for each axis) and return the space group name as used in GSAS
3616    '''
3617    rspc = spcgroup.replace(' ','').upper()
3618    # deal with rhombohedral and hexagonal setting designations
3619    rhomb =