source: trunk/GSASIIspc.py @ 3795

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

fix space group & incommensurate space group displays for mag. incomensurate cases
disable cell refine in Unit Cell List if no peaks
another shot at incommensurate mag. structure factor calc.

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