source: trunk/GSASIIspc.py @ 3961

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

set self.tutorialPath = '.\G2tutorial' if others fail. Don't know if this is best?
Add S301 to trigonal to 3m1 list of generalized mustrain terms - seems to improve things but plot is weird. Might need more work.
use get to fix missing Z problem
set up getting better map density from density slice instead of 3/4D full map

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