source: trunk/GSASIIspc.py @ 3754

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

begin work on incommensurate mag. structure factors
fix problem of magnetic reflection extinctions
replace analytical derivatives of magnetic moments with numerical ones using deltM = 0.0001

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