source: trunk/GSASIIspc.py

Last change on this file was 5171, checked in by vondreele, 8 months ago

fix ApplyModeDisp? to follow site symmetry rules.
Add xyz pointers to G2spc.CSXinel arrays for special positions
Cleanup effect in G2strIO because GetCSxinel now return 3 arrays.

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