source: trunk/GSASIIspc.py @ 3729

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

make super group ComboBox? allow editing for nonOSX.

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