source: trunk/GSASIIspc.py @ 4651

Last change on this file since 4651 was 4651, checked in by vondreele, 15 months ago

problem with "generalized" mustrain for R3 now fixed.The S301 term should not have been present.

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