source: trunk/GSASIIspc.py @ 4213

Last change on this file since 4213 was 4213, checked in by toby, 2 years ago

multiple small changes to allow docs build under 3.x

  • 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: 2019-12-20 05:04:03 +0000 (Fri, 20 Dec 2019) $
12# $Author: toby $
13# $Revision: 4213 $
14# $URL: trunk/GSASIIspc.py $
15# $Id: GSASIIspc.py 4213 2019-12-20 05:04:03Z toby $
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: 4213 $")
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)+10.)%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    cell0 = np.zeros(3,dtype=np.int32)
2086    if Move:
2087        X,cell0 = MoveToUnitCell(X)
2088    for ic,cen in enumerate(icen):
2089        C = np.array(cen)
2090        for invers in range(inv):
2091            for io,[M,T] in enumerate(SGData['SGOps']):
2092                idup = ((io+1)+100*ic)*(1-2*invers)
2093                XT = np.inner(M,X)+T
2094                if len(Uij):
2095                    U = Uij2U(Uij)
2096                    U = np.inner(M,np.inner(U,M).T)
2097                    newUij = U2Uij(U)
2098                if invers:
2099                    XT = -XT
2100                XT += C
2101                cell = np.zeros(3,dtype=np.int32)+cell0
2102                cellj = np.zeros(3,dtype=np.int32)
2103                if Move:
2104                    newX,cellj = MoveToUnitCell(XT)
2105                else:
2106                    newX = XT
2107                cell += cellj
2108                if All:
2109                    if np.allclose(newX,X,atol=0.0002):
2110                        idup = False
2111                else:
2112                    if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
2113                        idup = False
2114                if All or idup:
2115                    XYZEquiv.append(newX)
2116                    Idup.append(idup)
2117                    Cell.append(cell)
2118                    if len(Uij):
2119                        UijEquiv.append(newUij)
2120                    if len(SpnFlp):
2121                        spnflp.append(SpnFlp[mj])
2122                    else:
2123                        spnflp.append(1)
2124                mj += 1
2125    if len(Uij):
2126        return zip(XYZEquiv,UijEquiv,Idup,Cell,spnflp)
2127    else:
2128        return zip(XYZEquiv,Idup,Cell,spnflp)
2129       
2130def GenHKL(HKL,SGData):
2131    ''' Generates all equivlent reflections including Friedel pairs
2132    :param HKL:  [h,k,l] must be integral values
2133    :param SGData: space group data obtained from SpcGroup
2134    :returns: array Uniq: equivalent reflections
2135    '''
2136   
2137    Ops = SGData['SGOps']
2138    OpM = np.array([op[0] for op in Ops])
2139    Uniq = np.inner(OpM,HKL)
2140    Uniq = list(Uniq)+list(-1*Uniq)
2141    return np.array(Uniq)
2142
2143def GenHKLf(HKL,SGData):
2144    '''
2145    Uses old GSAS Fortran routine genhkl.for
2146
2147    :param HKL:  [h,k,l] must be integral values for genhkl.for to work
2148    :param SGData: space group data obtained from SpcGroup
2149    :returns: iabsnt,mulp,Uniq,phi
2150
2151     *   iabsnt = True if reflection is forbidden by symmetry
2152     *   mulp = reflection multiplicity including Friedel pairs
2153     *   Uniq = numpy array of equivalent hkl in descending order of h,k,l
2154     *   phi = phase offset for each equivalent h,k,l
2155
2156    '''
2157    hklf = list(HKL)+[0,]       #could be numpy array!
2158    Ops = SGData['SGOps']
2159    OpM = np.array([op[0] for op in Ops],order='F')
2160    OpT = np.array([op[1] for op in Ops])
2161    Cen = np.array([cen for cen in SGData['SGCen']],order='F')
2162   
2163    import pyspg
2164    Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
2165    h,k,l,f = Uniq
2166    Uniq=np.array(list(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq])))
2167    phi = f[:Nuniq]
2168    return iabsnt,mulp,Uniq,phi
2169   
2170def checkSSLaue(HKL,SGData,SSGData):
2171    #Laue check here - Toss HKL if outside unique Laue part
2172    h,k,l,m = HKL
2173    if SGData['SGLaue'] == '2/m':
2174        if SGData['SGUniq'] == 'a':
2175            if 'a' in SSGData['modSymb'] and h == 0 and m < 0:
2176                return False
2177            elif 'b' in SSGData['modSymb'] and k == 0 and l ==0 and m < 0:
2178                return False
2179            else:
2180                return True
2181        elif SGData['SGUniq'] == 'b':
2182            if 'b' in SSGData['modSymb'] and k == 0 and m < 0:
2183                return False
2184            elif 'a' in SSGData['modSymb'] and h == 0 and l ==0 and m < 0:
2185                return False
2186            else:
2187                return True
2188        elif SGData['SGUniq'] == 'c':
2189            if 'g' in SSGData['modSymb'] and l == 0 and m < 0:
2190                return False
2191            elif 'a' in SSGData['modSymb'] and h == 0 and k ==0 and m < 0:
2192                return False
2193            else:
2194                return True
2195    elif SGData['SGLaue'] == 'mmm':
2196        if 'a' in SSGData['modSymb']:
2197            if h == 0 and m < 0:
2198                return False
2199            else:
2200                return True
2201        elif 'b' in SSGData['modSymb']:
2202            if k == 0 and m < 0:
2203                return False
2204            else:
2205                return True
2206        elif 'g' in SSGData['modSymb']:
2207            if l == 0 and m < 0:
2208                return False
2209            else:
2210                return True
2211    else:   #tetragonal, trigonal, hexagonal (& triclinic?)
2212        if l == 0 and m < 0:
2213            return False
2214        else:
2215            return True
2216       
2217def checkHKLextc(HKL,SGData):
2218    '''
2219    Checks if reflection extinct - does not check centering
2220
2221    :param HKL:  [h,k,l]
2222    :param SGData: space group data obtained from SpcGroup
2223    :returns: True if extinct; False if allowed
2224
2225    '''
2226    Ops = SGData['SGOps']
2227    OpM = np.array([op[0] for op in Ops])
2228    OpT = np.array([op[1] for op in Ops])
2229    HKLS = np.array([HKL,-HKL])     #Freidel's Law
2230    DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,3))
2231    PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2232    for dhkl,phkl in zip(DHKL,PHKL)[1:]:    #skip identity
2233        if dhkl.any():
2234            continue
2235        else:
2236            if phkl%1.:
2237                return True
2238    return False
2239
2240def checkMagextc(HKL,SGData):
2241    '''
2242    Checks if reflection magnetically extinct; does fullcheck (centering, too)
2243    uses algorthm from Gallego, et al., J. Appl. Cryst. 45, 1236-1247 (2012)
2244
2245    :param HKL:  [h,k,l]
2246    :param SGData: space group data obtained from SpcGroup; must have magnetic symmetry SpnFlp data
2247    :returns: True if magnetically extinct; False if allowed (to match GenHKLf)
2248
2249    '''
2250    Ops = SGData['SGOps']
2251    Ncen = len(SGData['SGCen'])
2252    OpM = np.array([op[0] for op in Ops])
2253    OpT = np.array([op[1] for op in Ops])
2254    if SGData['SGInv'] and not SGData['SGFixed']:
2255        OpM = np.vstack((OpM,-OpM))
2256        OpT = np.vstack((OpT,-OpT))%1.
2257    OpM = np.reshape(np.array(list(OpM)*Ncen),(-1,3,3))
2258    OpT = np.reshape(np.array([OpT+cen for cen in SGData['SGCen']]),(-1,3))
2259    Spn = SGData['SpnFlp'][:len(OpM)]
2260    Mag = np.array([nl.det(opm) for opm in OpM])*Spn
2261    DHKL = np.reshape(np.inner(HKL,OpM),(-1,3))
2262    PHKL = np.reshape(np.cos(2.0*np.pi*np.inner(HKL,OpT))*Mag,(-1,))[:,nxs,nxs]*OpM     #compute T(R,theta) eq(7)
2263    Ftest = np.random.rand(3)       #random magnetic moment
2264    Psum = np.zeros(3)
2265    nsum = 0.
2266    nA = 0
2267    for dhkl,phkl in zip(DHKL,PHKL):
2268        if not np.allclose(dhkl,HKL):           #test for eq(5)
2269            continue
2270        else:
2271            nA += 1
2272            nsum += np.trace(phkl)          #eq(8)
2273            pterm = np.inner(Ftest,phkl)    #eq(9)
2274            Psum += pterm
2275    if nsum/nA > 1.:        #only need to look at nA=1 frok eq(8)
2276        return False
2277    if np.allclose(Psum,np.zeros(3)):
2278        return True
2279    else:
2280        if np.inner(HKL,Psum):
2281            return True
2282        return False
2283   
2284def checkSSextc(HKL,SSGData):
2285    Ops = SSGData['SSGOps']
2286    OpM = np.array([op[0] for op in Ops])
2287    OpT = np.array([op[1] for op in Ops])
2288    HKLS = np.array([HKL,-HKL])     #Freidel's Law
2289    DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,4))
2290    PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2291    for dhkl,phkl in list(zip(DHKL,PHKL))[1:]:    #skip identity
2292        if dhkl.any():
2293            continue
2294        else:
2295            if phkl%1.:
2296                return False
2297    return True
2298   
2299################################################################################
2300#### Site symmetry tables
2301################################################################################
2302     
2303OprName = {
2304    '-6643':       ['-1',1],'6479' :    ['2(z)',2],'-6479':     ['m(z)',3],
2305    '6481' :     ['m(y)',4],'-6481':    ['2(y)',5],'6641' :     ['m(x)',6],
2306    '-6641':     ['2(x)',7],'6591' :  ['m(+-0)',8],'-6591':   ['2(+-0)',9],
2307    '6531' :  ['m(110)',10],'-6531': ['2(110)',11],'6537' :    ['4(z)',12],
2308    '-6537':   ['-4(z)',13],'975'  : ['3(111)',14],'6456' :       ['3',15],
2309    '-489' :  ['3(+--)',16],'483'  : ['3(-+-)',17],'-969' :  ['3(--+)',18],
2310    '819'  :  ['m(+0-)',19],'-819' : ['2(+0-)',20],'2431' :  ['m(0+-)',21],
2311    '-2431':  ['2(0+-)',22],'-657' :  ['m(xz)',23],'657'  :   ['2(xz)',24],
2312    '1943' :   ['-4(x)',25],'-1943':   ['4(x)',26],'-2429':   ['m(yz)',27],
2313    '2429' :   ['2(yz)',28],'639'  :  ['-4(y)',29],'-639' :    ['4(y)',30],
2314    '-6484':   ['2(010)',4],'6484' :  ['m(010)',5],'-6668':   ['2(100)',6],
2315    '6668' :   ['m(100)',7],'-6454': ['2(120)',18],'6454' :  ['m(120)',19],
2316    '-6638':  ['2(210)',20],'6638' : ['m(210)',21],   #search in SytSym ends at m(210)
2317    '2223' : ['3(+++)2',39],
2318    '6538' :   ['6(z)1',40],'-2169':['3(--+)2',41],'2151' : ['3(+--)2',42],
2319    '2205' :['-3(-+-)2',43],'-2205':[' (-+-)2',44],'489'  :['-3(+--)1',45],
2320    '801'  :   ['4(y)1',46],'1945' :  ['4(x)3',47],'-6585': ['-4(z)3 ',48],
2321    '6585' :   ['4(z)3',49],'6584' :  ['3(z)2',50],'6666' :  ['6(z)5 ',51],
2322    '6643' :       ['1',52],'-801' : ['-4(y)1',53],'-1945': ['-4(x)3 ',54],
2323    '-6666':  ['-6(z)5',55],'-6538': ['-6(z)1',56],'-2223':['-3(+++)2',57],
2324    '-975' :['-3(+++)1',58],'-6456': ['-3(z)1',59],'-483' :['-3(-+-)1',60],
2325    '969'  :['-3(--+)1',61],'-6584': ['-3(z)2',62],'2169' :['-3(--+)2',63],
2326    '-2151':['-3(+--)2',64],   }                               
2327
2328KNsym = {
2329    '0'         :'    1   ','1'         :'   -1   ','64'        :'    2(x)','32'        :'    m(x)',
2330    '97'        :'  2/m(x)','16'        :'    2(y)','8'         :'    m(y)','25'        :'  2/m(y)',
2331    '2'         :'    2(z)','4'         :'    m(z)','7'         :'  2/m(z)','134217728' :'   2(yz)',
2332    '67108864'  :'   m(yz)','201326593' :' 2/m(yz)','2097152'   :'  2(0+-)','1048576'   :'  m(0+-)',
2333    '3145729'   :'2/m(0+-)','8388608'   :'   2(xz)','4194304'   :'   m(xz)','12582913'  :' 2/m(xz)',
2334    '524288'    :'  2(+0-)','262144'    :'  m(+0-)','796433'    :'2/m(+0-)','1024'      :'   2(xy)',
2335    '512'       :'   m(xy)','1537'      :' 2/m(xy)','256'       :'  2(+-0)','128'       :'  m(+-0)',
2336    '385'       :'2/m(+-0)','76'        :'  mm2(x)','52'        :'  mm2(y)','42'        :'  mm2(z)',
2337    '135266336' :' mm2(yz)','69206048'  :'mm2(0+-)','8650760'   :' mm2(xz)','4718600'   :'mm2(+0-)',
2338    '1156'      :' mm2(xy)','772'       :'mm2(+-0)','82'        :'  222   ','136314944' :'  222(x)',
2339    '8912912'   :'  222(y)','1282'      :'  222(z)','127'       :'  mmm   ','204472417' :'  mmm(x)',
2340    '13369369'  :'  mmm(y)','1927'      :'  mmm(z)','33554496'  :'  4(x)','16777280'  :' -4(x)',
2341    '50331745'  :'4/m(x)'  ,'169869394' :'422(x)','84934738'  :'-42m(x)','101711948' :'4mm(x)',
2342    '254804095' :'4/mmm(x)','536870928 ':'  4(y)','268435472' :' -4(y)','805306393' :'4/m(y)',
2343    '545783890' :'422(y)','272891986' :'-42m(y)','541327412' :'4mm(y)','818675839' :'4/mmm(y)',
2344    '2050'      :'  4(z)','4098'      :' -4(z)','6151'      :'4/m(z)','3410'      :'422(z)',
2345    '4818'      :'-42m(z)','2730'      :'4mm(z)','8191'      :'4/mmm(z)','8192'      :'  3(111)',
2346    '8193'      :' -3(111)','2629888'   :' 32(111)','1319040'   :' 3m(111)','3940737'   :'-3m(111)',
2347    '32768'     :'  3(+--)','32769'     :' -3(+--)','10519552'  :' 32(+--)','5276160'   :' 3m(+--)',
2348    '15762945'  :'-3m(+--)','65536'     :'  3(-+-)','65537'     :' -3(-+-)','134808576' :' 32(-+-)',
2349    '67437056'  :' 3m(-+-)','202180097' :'-3m(-+-)','131072'    :'  3(--+)','131073'    :' -3(--+)',
2350    '142737664' :' 32(--+)','71434368'  :' 3m(--+)','214040961' :'-3m(--+)','237650'    :'   23   ',
2351    '237695'    :'   m3   ','715894098' :'   432  ','358068946' :'  -43m  ','1073725439':'   m3m  ',
2352    '68157504'  :' mm2(d100)','4456464'   :' mm2(d010)','642'       :' mm2(d001)','153092172' :'-4m2(x)',
2353    '277348404' :'-4m2(y)','5418'      :'-4m2(z)','1075726335':'  6/mmm ','1074414420':'-6m2(100)',
2354    '1075070124':'-6m2(120)','1075069650':'   6mm  ','1074414890':'   622  ','1073758215':'   6/m  ',
2355    '1073758212':'   -6   ','1073758210':'    6   ','1073759865':'-3m(100)','1075724673':'-3m(120)',
2356    '1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
2357    '1073758209':'   -3   ','1073758208':'    3   ','1074135143':'mmm(100)','1075314719':'mmm(010)',
2358    '1073743751':'mmm(110)','1074004034':' mm2(z100)','1074790418':' mm2(z010)','1073742466':' mm2(z110)',
2359    '1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
2360    '1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
2361    '1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
2362    '1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
2363    '1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
2364    '1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
2365    '1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
2366    '1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':'   -1   ',
2367    }
2368
2369NXUPQsym = {
2370    '1'        :(28,29,28,28),'-1'       :( 1,29,28, 0),'2(x)'     :(12,18,12,25),'m(x)'     :(25,18,12,25),
2371    '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),
2372    '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),
2373    '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),
2374    '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),
2375    '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),
2376    '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),
2377    '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),
2378    'mm2(yz)'  :(10,13, 0,-1),'mm2(0+-)' :(11,13, 0,-1),'mm2(xz)'  :( 8,12, 0,-1),'mm2(+0-)' :( 9,12, 0,-1),
2379    'mm2(xy)'  :( 6,11, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222'      :( 1,10, 0,-1),'222(x)'   :( 1,13, 0,-1),
2380    '222(y)'   :( 1,12, 0,-1),'222(z)'   :( 1,11, 0,-1),'mmm'      :( 1,10, 0,-1),'mmm(x)'   :( 1,13, 0,-1),
2381    'mmm(y)'   :( 1,12, 0,-1),'mmm(z)'   :( 1,11, 0,-1),'4(x)'     :(12, 4,12, 0),'-4(x)'    :( 1, 4,12, 0),
2382    '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),
2383    '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),
2384    '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,),
2385    '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),
2386    '-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),
2387    '-3(111)'  :( 1, 5, 2, 0),'32(111)'  :( 1, 5, 0, 2),'3m(111)'  :( 2, 5, 0, 2),'-3m(111)' :( 1, 5, 0,-1),
2388    '3(+--)'   :( 5, 8, 5, 0),'-3(+--)'  :( 1, 8, 5, 0),'32(+--)'  :( 1, 8, 0, 5),'3m(+--)'  :( 5, 8, 0, 5),
2389    '-3m(+--)' :( 1, 8, 0,-1),'3(-+-)'   :( 4, 7, 4, 0),'-3(-+-)'  :( 1, 7, 4, 0),'32(-+-)'  :( 1, 7, 0, 4),
2390    '3m(-+-)'  :( 4, 7, 0, 4),'-3m(-+-)' :( 1, 7, 0,-1),'3(--+)'   :( 3, 6, 3, 0),'-3(--+)'  :( 1, 6, 3, 0),
2391    '32(--+)'  :( 1, 6, 0, 3),'3m(--+)'  :( 3, 6, 0, 3),'-3m(--+)' :( 1, 6, 0,-1),'23'       :( 1, 1, 0, 0),
2392    'm3'       :( 1, 1, 0, 0),'432'      :( 1, 1, 0, 0),'-43m'     :( 1, 1, 0, 0),'m3m'      :( 1, 1, 0, 0),
2393    'mm2(d100)':(12,13, 0,-1),'mm2(d010)':(13,12, 0,-1),'mm2(d001)':(14,11, 0,-1),'-4m2(x)'  :( 1, 4, 0,-1),
2394    '-4m2(y)'  :( 1, 3, 0,-1),'-4m2(z)'  :( 1, 2, 0,-1),'6/mmm'    :( 1, 9, 0,-1),'-6m2(100)':( 1, 9, 0,-1),
2395    '-6m2(120)':( 1, 9, 0,-1),'6mm'      :(14, 9, 0,-1),'622'      :( 1, 9, 0,-1),'6/m'      :( 1, 9,14,-1),
2396    '-6'       :( 1, 9,14, 0),'6'        :(14, 9,14, 0),'-3m(100)' :( 1, 9, 0,-1),'-3m(120)' :( 1, 9, 0,-1),
2397    '3m(100)'  :(14, 9, 0,14),'3m(120)'  :(14, 9, 0,14),'32(100)'  :( 1, 9, 0,14),'32(120)'  :( 1, 9, 0,14),
2398    '-3'       :( 1, 9,14, 0),'3'        :(14, 9,14, 0),'mmm(100)' :( 1,14, 0,-1),'mmm(010)' :( 1,15, 0,-1),
2399    'mmm(110)' :( 1,11, 0,-1),'mm2(z100)':(14,14, 0,-1),'mm2(z010)':(14,15, 0,-1),'mm2(z110)':(14,11, 0,-1),
2400    'mm2(100)' :(12,14, 0,-1),'mm2(010)' :(13,15, 0,-1),'mm2(110)' :( 6,11, 0,-1),'mm2(120)' :(15,14, 0,-1),
2401    'mm2(210)' :(16,15, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222(100)' :( 1,14, 0,-1),'222(010)' :( 1,15, 0,-1),
2402    '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),
2403    '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),
2404    'm(001)'   :(23,16,14,23),'m(100)'   :(26,25,12,26),'m(010)'   :(27,28,13,27),'m(110)'   :(18,19, 6,18),
2405    'm(120)'   :(24,27,15,24),'m(210)'   :(25,26,16,25),'m(+-0)'   :(17,20, 7,17),'2(001)'   :(14,16,14,23),
2406    '2(100)'   :(12,25,12,26),'2(010)'   :(13,28,13,27),'2(110)'   :( 6,19, 6,18),'2(120)'   :(15,27,15,24),
2407    '2(210)'   :(16,26,16,25),'2(+-0)'   :( 7,20, 7,17),'-1'       :( 1,29,28, 0)
2408    }
2409       
2410CSxinel = [[],      # 0th empty - indices are Fortran style
2411    [[0,0,0],[ 0.0, 0.0, 0.0]],      #1  0  0  0
2412    [[1,1,1],[ 1.0, 1.0, 1.0]],      #2  X  X  X
2413    [[1,1,1],[ 1.0, 1.0,-1.0]],      #3  X  X -X
2414    [[1,1,1],[ 1.0,-1.0, 1.0]],      #4  X -X  X
2415    [[1,1,1],[ 1.0,-1.0,-1.0]],      #5 -X  X  X
2416    [[1,1,0],[ 1.0, 1.0, 0.0]],      #6  X  X  0
2417    [[1,1,0],[ 1.0,-1.0, 0.0]],      #7  X -X  0
2418    [[1,0,1],[ 1.0, 0.0, 1.0]],      #8  X  0  X
2419    [[1,0,1],[ 1.0, 0.0,-1.0]],      #9  X  0 -X
2420    [[0,1,1],[ 0.0, 1.0, 1.0]],      #10  0  Y  Y
2421    [[0,1,1],[ 0.0, 1.0,-1.0]],      #11 0  Y -Y
2422    [[1,0,0],[ 1.0, 0.0, 0.0]],      #12  X  0  0
2423    [[0,1,0],[ 0.0, 1.0, 0.0]],      #13  0  Y  0
2424    [[0,0,1],[ 0.0, 0.0, 1.0]],      #14  0  0  Z
2425    [[1,1,0],[ 1.0, 2.0, 0.0]],      #15  X 2X  0
2426    [[1,1,0],[ 2.0, 1.0, 0.0]],      #16 2X  X  0
2427    [[1,1,2],[ 1.0, 1.0, 1.0]],      #17  X  X  Z
2428    [[1,1,2],[ 1.0,-1.0, 1.0]],      #18  X -X  Z
2429    [[1,2,1],[ 1.0, 1.0, 1.0]],      #19  X  Y  X
2430    [[1,2,1],[ 1.0, 1.0,-1.0]],      #20  X  Y -X
2431    [[1,2,2],[ 1.0, 1.0, 1.0]],      #21  X  Y  Y
2432    [[1,2,2],[ 1.0, 1.0,-1.0]],      #22  X  Y -Y
2433    [[1,2,0],[ 1.0, 1.0, 0.0]],      #23  X  Y  0
2434    [[1,0,2],[ 1.0, 0.0, 1.0]],      #24  X  0  Z
2435    [[0,1,2],[ 0.0, 1.0, 1.0]],      #25  0  Y  Z
2436    [[1,1,2],[ 1.0, 2.0, 1.0]],      #26  X 2X  Z
2437    [[1,1,2],[ 2.0, 1.0, 1.0]],      #27 2X  X  Z
2438    [[1,2,3],[ 1.0, 1.0, 1.0]],      #28  X  Y  Z
2439    ]
2440
2441CSuinel = [[],      # 0th empty - indices are Fortran style
2442    [[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
2443    [[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
2444    [[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
2445    [[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
2446    [[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
2447    [[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
2448    [[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
2449    [[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
2450    [[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
2451    [[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
2452    [[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
2453    [[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
2454    [[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
2455    [[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
2456    [[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
2457    [[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
2458    [[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
2459    [[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
2460    [[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
2461    [[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
2462    [[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
2463    [[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
2464    [[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
2465    [[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
2466    [[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
2467    [[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
2468    [[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
2469    [[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
2470    [[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
2471    ]
2472   
2473################################################################################
2474#### Site symmetry routines
2475################################################################################
2476   
2477def GetOprPtrName(key):
2478    'Needs a doc string'
2479    try:
2480        oprName = OprName[key][0]
2481    except KeyError:
2482        return key
2483    return oprName.replace('(','').replace(')','')
2484
2485def GetOprPtrNumber(key):
2486    'Needs a doc string'
2487    try:
2488        return OprName[key][1]
2489    except KeyError:
2490        return key
2491
2492def GetOprName(key):
2493    'Needs a doc string'
2494    return OprName[key][0]
2495
2496def GetKNsym(key):
2497    'Needs a doc string'
2498    try:
2499        return KNsym[key].strip()
2500    except KeyError:
2501        return 'sp'
2502
2503def GetNXUPQsym(siteSym):
2504    '''       
2505    The codes XUPQ are for lookup of symmetry constraints for position(X), thermal parm(U) & magnetic moments (P & Q)
2506    '''
2507    return NXUPQsym[siteSym]
2508
2509def GetCSxinel(siteSym): 
2510    "returns Xyz terms, multipliers, GUI flags"
2511    indx = GetNXUPQsym(siteSym.strip())
2512    return CSxinel[indx[0]]
2513   
2514def GetCSuinel(siteSym):
2515    "returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
2516    indx = GetNXUPQsym(siteSym.strip())
2517    return CSuinel[indx[1]]
2518   
2519def GetCSpqinel(SpnFlp,dupDir): 
2520    "returns Mxyz terms, multipliers, GUI flags"
2521    CSI = [[1,2,3],[1.0,1.0,1.0]]
2522    for sopr in dupDir:
2523#        print (sopr,dupDir[sopr])
2524        opr = sopr.replace("'",'')
2525        indx = GetNXUPQsym(opr)
2526        if SpnFlp[dupDir[sopr]] > 0:
2527            csi = CSxinel[indx[2]]  #P
2528        else:
2529            csi = CSxinel[indx[3]]  #Q
2530#        print(opr,indx,csi,CSI)
2531        if not len(csi):
2532            return [[0,0,0],[0.,0.,0.]]
2533        for kcs in [0,1,2]:
2534            if csi[0][kcs] == 0 and CSI[0][kcs] != 0:
2535                jcs = CSI[0][kcs]
2536                for ics in [0,1,2]:
2537                    if CSI[0][ics] == jcs:
2538                        CSI[0][ics] = 0
2539                        CSI[1][ics] = 0.
2540                    elif CSI[0][ics] > jcs:
2541                        CSI[0][ics] = CSI[0][ics]-1
2542            elif (CSI[0][kcs] == csi[0][kcs]) and (CSI[1][kcs] != csi[1][kcs]):
2543                CSI[1][kcs] = csi[1][kcs]
2544            elif CSI[0][kcs] >= csi[0][kcs]:
2545                CSI[0][kcs] = min(CSI[0][kcs],csi[0][kcs])
2546                if CSI[1][kcs] != csi[1][kcs]:
2547                    if CSI[1][kcs] == 1.:
2548                        CSI[1][kcs] = csi[1][kcs]
2549#        print(CSI)
2550    return CSI
2551   
2552def getTauT(tau,sop,ssop,XYZ,wave=np.zeros(3)):
2553    phase = np.sum(XYZ*wave)
2554    ssopinv = nl.inv(ssop[0])
2555    mst = ssopinv[3][:3]
2556    epsinv = ssopinv[3][3]
2557    sdet = nl.det(sop[0])
2558    ssdet = nl.det(ssop[0])
2559    dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
2560    dT = 1.0
2561    if np.any(dtau%.5):
2562        sumdtau = np.sum(dtau%.5)
2563        dT = 0.
2564        if np.abs(sumdtau-.5) > 1.e-4:
2565            dT = np.tan(np.pi*sumdtau)
2566    tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]+phase)
2567    return sdet,ssdet,dtau,dT,tauT
2568   
2569def OpsfromStringOps(A,SGData,SSGData):
2570    SGOps = SGData['SGOps']
2571    SSGOps = SSGData['SSGOps']
2572    Ax = A.split('+')
2573    Ax[0] = int(Ax[0])
2574    iC = 1
2575    if Ax[0] < 0:
2576        iC = -1
2577    iAx = abs(Ax[0])
2578    nA = iAx%100-1
2579    nC = iAx//100
2580    unit = [0,0,0]
2581    if len(Ax) > 1:
2582        unit = eval('['+Ax[1]+']')
2583    return SGOps[nA],SSGOps[nA],iC,SGData['SGCen'][nC],unit
2584   
2585def GetSSfxuinel(waveType,Stype,nH,XYZ,SGData,SSGData,debug=False):
2586   
2587    def orderParms(CSI):
2588        parms = [0,]
2589        for csi in CSI:
2590            for i in [0,1,2]:
2591                if csi[i] not in parms:
2592                    parms.append(csi[i])
2593        for csi in CSI:
2594            for i in [0,1,2]:
2595                csi[i] = parms.index(csi[i])
2596        return CSI
2597       
2598    def fracCrenel(tau,Toff,Twid):
2599        Tau = (tau-Toff[:,nxs])%1.
2600        A = np.where(Tau<Twid[:,nxs],1.,0.)
2601        return A
2602       
2603    def fracFourier(tau,nH,fsin,fcos):
2604        SA = np.sin(2.*nH*np.pi*tau)
2605        CB = np.cos(2.*nH*np.pi*tau)
2606        A = SA[nxs,nxs,:]*fsin[:,:,nxs]
2607        B = CB[nxs,nxs,:]*fcos[:,:,nxs]
2608        return A+B
2609       
2610    def posFourier(tau,nH,psin,pcos):
2611        SA = np.sin(2*nH*np.pi*tau)
2612        CB = np.cos(2*nH*np.pi*tau)
2613        A = SA[nxs,nxs,:]*psin[:,:,nxs]
2614        B = CB[nxs,nxs,:]*pcos[:,:,nxs]
2615        return A+B   
2616
2617    def posZigZag(tau,Tmm,XYZmax):
2618        DT = Tmm[1]-Tmm[0]
2619        slopeUp = 2.*XYZmax/DT
2620        slopeDn = 2.*XYZmax/(1.-DT)
2621        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])
2622        return A
2623
2624    def posBlock(tau,Tmm,XYZmax):
2625        A = np.array([np.where(Tmm[0] < t <= Tmm[1],XYZmax,-XYZmax) for t in tau])
2626        return A
2627       
2628    def DoFrac():
2629        delt2 = np.eye(2)*0.001
2630        dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2631        dFTP = []
2632        if siteSym == '1':
2633            CSI = [[1,0],[2,0]],2*[[1.,0.],]
2634        elif siteSym == '-1':
2635            CSI = [[1,0],[0,0]],2*[[1.,0.],]
2636        else:
2637            FSC = np.ones(2,dtype='i')
2638            CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
2639            if 'Crenel' in waveType:
2640                dF = np.zeros_like(tau)
2641            else:
2642                dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2643            dFT = np.zeros_like(dF)
2644            dFTP = []
2645            for i in SdIndx:
2646                sop = Sop[i]
2647                ssop = SSop[i]           
2648                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2649                fsc = np.ones(2,dtype='i')
2650                if 'Crenel' in waveType:
2651                    dFT = np.zeros_like(tau)
2652                    fsc = [1,1]
2653                else:   #Fourier
2654                    dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
2655                    dFT = nl.det(sop[0])*dFT
2656                    dFT = dFT[:,np.argsort(tauT)]
2657                    dFT[0] *= ssdet
2658                    dFT[1] *= sdet
2659                    dFTP.append(dFT)
2660               
2661                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2662                        fsc = [1,1]
2663                        if dT:
2664                            CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
2665                        else:
2666                            CSI = [[[1,0],[0,0]],[[1.,0.],[0.,0.]]]
2667                        FSC = np.zeros(2,dtype='i')
2668                        return CSI,dF,dFTP
2669                    else:
2670                        for i in range(2):
2671                            if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
2672                                fsc[i] = 1
2673                            else:
2674                                fsc[i] = 0
2675                        FSC &= fsc
2676                        if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc)
2677            n = -1
2678            for i,F in enumerate(FSC):
2679                if F:
2680                    n += 1
2681                    CSI[0][i] = n+1
2682                    CSI[1][i] = 1.0
2683           
2684        return CSI,dF,dFTP
2685       
2686    def DoXYZ():
2687        delt5 = np.ones(5)*0.001
2688        delt6 = np.eye(6)*0.001
2689        if 'Fourier' in waveType:
2690            dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2691              #3x6x12 modulated position array (X,Spos,tau)& force positive
2692        elif waveType in ['ZigZag','Block']:
2693            if waveType == 'ZigZag':
2694                dX = posZigZag(tau,delt5[:2],delt5[2:])
2695            else:
2696                dX = posBlock(tau,delt5[:2],delt5[2:])
2697        dXTP = []
2698        if siteSym == '1':
2699            CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2700        elif siteSym == '-1':
2701            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.],]
2702        else:
2703            if 'Fourier' in waveType:
2704                CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2705            elif waveType in ['ZigZag','Block']:
2706                CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2707                    np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2708            XSC = np.ones(6,dtype='i')
2709            dXTP = []
2710            for i in SdIndx:
2711                sop = Sop[i]
2712                ssop = SSop[i]
2713                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2714                xsc = np.ones(6,dtype='i')
2715                if 'Fourier' in waveType:
2716                    dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,nxs,nxs]
2717                elif waveType == 'ZigZag':
2718                    dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2719                elif waveType == 'Block':
2720                    dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2721                dXT = np.inner(sop[0],dXT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2722                dXT = np.swapaxes(dXT,1,2)      # back to array(3x6x49)
2723                dXT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2724                dXTP.append(dXT)
2725                if waveType == 'Fourier':
2726                    for i in range(3):
2727                        if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2728                            xsc[i] = 0
2729                        if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2730                            xsc[i+3] = 0
2731                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2732                        xsc[3:6] = 0
2733                        CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2734                            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2735                        if dT:
2736                            if '(x)' in siteSym:
2737                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2738                                if 'm' in siteSym and len(SdIndx) == 1:
2739                                    CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2740                            elif '(y)' in siteSym:
2741                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2742                                if 'm' in siteSym and len(SdIndx) == 1:
2743                                    CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2744                            elif '(z)' in siteSym:
2745                                CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2746                                if 'm' in siteSym and len(SdIndx) == 1:
2747                                    CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2748                        else:
2749                            CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2750                    if '4/mmm' in laue:
2751                        if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2752                            if '(xy)' in siteSym:
2753                                CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2754                                if dT:
2755                                    CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2756                                else:
2757                                    CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2758                        if '(xy)' in siteSym or '(+-0)' in siteSym:
2759                            mul = 1
2760                            if '(+-0)' in siteSym:
2761                                mul = -1
2762                            if np.allclose(dX[0,0,:],dXT[1,0,:]):
2763                                CSI[0][3:5] = [[11,0,0],[11,0,0]]
2764                                CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2765                                xsc[3:5] = 0
2766                            if np.allclose(dX[0,3,:],dXT[0,4,:]):
2767                                CSI[0][:2] = [[12,0,0],[12,0,0]]
2768                                CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2769                                xsc[:2] = 0
2770                else:
2771                    for i in range(3):
2772                        if not np.allclose(dX[:,i],dXT[i,:,i]):
2773                            xsc[i] = 0
2774                XSC &= xsc
2775                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2776            if waveType == 'Fourier':
2777                n = -1
2778                if debug: print (XSC)
2779                for i,X in enumerate(XSC):
2780                    if X:
2781                        n += 1
2782                        CSI[0][i][0] = n+1
2783                        CSI[1][i][0] = 1.0
2784           
2785        return list(CSI),dX,dXTP
2786       
2787    def DoUij():
2788        delt12 = np.eye(12)*0.0001
2789        dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
2790        dUTP = []
2791        if siteSym == '1':
2792            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2793                [7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2794        elif siteSym == '-1':
2795            CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],   \
2796                6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2797        else:
2798            CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2799            USC = np.ones(12,dtype='i')
2800            dUTP = []
2801            dtau = 0.
2802            for i in SdIndx:
2803                sop = Sop[i]
2804                ssop = SSop[i]
2805                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2806                usc = np.ones(12,dtype='i')
2807                dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
2808                dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3)    #convert dUT to 12x49x3x3
2809                dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2810                dUT = np.array(U2Uij(dUijT))    #convert to 6x12x49
2811                dUT = dUT[:,:,np.argsort(tauT)]
2812                dUT[:,:6,:] *=(ssdet*sdet)
2813                dUTP.append(dUT)
2814                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2815                    if dT:
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                        [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2820                    else:
2821                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2822                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2823                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2824                        [0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2825                    if 'mm2(x)' in siteSym and dT:
2826                        CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2827                        USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2828                    elif '(xy)' in siteSym and dT:
2829                        CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2830                            [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2831                        CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2832                        USC = [1,1,1,1,1,1,1,1,1,1,1,1]                             
2833                    elif '(x)' in siteSym and dT:
2834                        CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2835                    elif '(y)' in siteSym and dT:
2836                        CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2837                    elif '(z)' in siteSym and dT:
2838                        CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2839                    for i in range(6):
2840                        if not USC[i]:
2841                            CSI[0][i] = [0,0,0]
2842                            CSI[1][i] = [0.,0.,0.]
2843                            CSI[0][i+6] = [0,0,0]
2844                            CSI[1][i+6] = [0.,0.,0.]
2845                else:                       
2846                    for i in range(6):
2847                        if not np.allclose(dU[i,i,:],dUT[i,i,:]):  #sin part
2848                            usc[i] = 0
2849                        if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]):   #cos part
2850                            usc[i+6] = 0
2851                    if np.any(dUT[1,0,:]):
2852                        if '4/m' in siteSym:
2853                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2854                            if ssop[1][3]:
2855                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2856                                usc[9] = 1
2857                            else:
2858                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2859                                usc[9] = 0
2860                        elif '4' in siteSym:
2861                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2862                            CSI[0][:2] = [[11,0,0],[11,0,0]]
2863                            if ssop[1][3]:
2864                                CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2865                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2866                                usc[2] = 0
2867                                usc[8] = 0
2868                                usc[3] = 1
2869                                usc[9] = 1
2870                            else:
2871                                CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2872                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2873                                usc[2] = 1
2874                                usc[8] = 1
2875                                usc[3] = 0               
2876                                usc[9] = 0
2877                        elif 'xy' in siteSym or '+-0' in siteSym:
2878                            if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2879                                CSI[0][4:6] = [[12,0,0],[12,0,0]]
2880                                CSI[0][6:8] = [[11,0,0],[11,0,0]]
2881                                CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2882                                CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2883                                usc[4:6] = 0
2884                                usc[6:8] = 0
2885                           
2886                    if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2887                USC &= usc
2888            if debug: print (USC)
2889            if not np.any(dtau%.5):
2890                n = -1
2891                for i,U in enumerate(USC):
2892                    if U:
2893                        n += 1
2894                        CSI[0][i][0] = n+1
2895                        CSI[1][i][0] = 1.0
2896   
2897        return list(CSI),dU,dUTP
2898   
2899    def DoMag():
2900        delt6 = np.eye(6)*0.001
2901        dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,nxs,nxs]
2902        dMTP = []
2903        CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2904        if siteSym == '1':
2905            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2906        elif siteSym in ['-1','mmm',]:
2907            CSI = 3*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3*[[0.,0.,0.],]+3*[[1.,0.,0.],]
2908        elif siteSym in ['4(z)','422(z)']:
2909            CSI[0][0][0] = CSI[0][4][1] = 1
2910            CSI[1][0][0] = 1.0
2911            CSI[1][4][1] = -1.0
2912        elif siteSym in ['-4m2(z)','422(z)',]:
2913            CSI[0][5][0] = 1
2914            CSI[1][5][0] = 1.0
2915        elif siteSym in ['-32(100)','-3',]:
2916            CSI[0][2][0] = 1
2917            CSI[1][2][0] = 1.0
2918        elif siteSym in ['3',]:
2919            CSI[0][0][0] = CSI[0][3][0] = CSI[0][4][0] = 1
2920            CSI[1][0][0] = -np.sqrt(3.0)
2921            CSI[1][3][0] = 2.0
2922            CSI[1][4][0] = 1.0
2923        elif siteSym in ['622','2(100)','32(100)',]:
2924            CSI[0][0][0] = CSI[0][1][0] = CSI[0][3][0] = 1
2925            CSI[1][0][0] = 1.0
2926            CSI[1][1][0] = 2.0
2927            CSI[1][3][0] = np.sqrt(3.0)
2928        else:
2929              #3x6x12 modulated moment array (M,Spos,tau)& force positive
2930            CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2931            MSC = np.ones(6,dtype='i')
2932            dMTP = []
2933            for i in SdIndx:
2934                sop = Sop[i]
2935                ssop = SSop[i]
2936                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2937                msc = np.ones(6,dtype='i')
2938                dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,nxs,nxs]
2939                dMT = np.inner(sop[0],dMT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2940                dMT = np.swapaxes(dMT,1,2)      # back to array(3x6x49)
2941                dMT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2942                dMTP.append(dMT)
2943                for i in range(3):
2944                    if not np.allclose(dM[i,i,:],sdet*dMT[i,i,:]):
2945                        msc[i] = 0
2946                    if not np.allclose(dM[i,i+3,:],sdet*dMT[i,i+3,:]):
2947                        msc[i+3] = 0
2948                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2949                    msc[3:6] = 0
2950                    CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2951                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2952                    if dT:
2953                        if '(x)' in siteSym:
2954                            CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2955                            if 'm' in siteSym and len(SdIndx) == 1:
2956                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2957                        elif '(y)' in siteSym:
2958                            CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2959                            if 'm' in siteSym and len(SdIndx) == 1:
2960                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2961                        elif '(z)' in siteSym:
2962                            CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2963                            if 'm' in siteSym and len(SdIndx) == 1:
2964                                CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2965                    else:
2966                        CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2967                if '4/mmm' in laue:
2968                    if siteSym in ['4/mmm(z)',]:
2969                        CSI = 3*[[0,0,0],]+[[0,0,0],[0,0,0],[1,0,0]],3*[[0.,0.,0.],]+3*[[1.,0.,0.],]
2970                    if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2971                        if '(xy)' in siteSym:
2972                            CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2973                            if dT:
2974                                CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2975                            else:
2976                                CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2977                    if '(xy)' in siteSym or '(+-0)' in siteSym:
2978                        mul = 1
2979                        if '(+-0)' in siteSym:
2980                            mul = -1
2981                        if np.allclose(dM[0,0,:],dMT[1,0,:]):
2982                            CSI[0][3:5] = [[11,0,0],[11,0,0]]
2983                            CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2984                            msc[3:5] = 0
2985                        if np.allclose(dM[0,3,:],dMT[0,4,:]):
2986                            CSI[0][:2] = [[12,0,0],[12,0,0]]
2987                            CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2988                            msc[:2] = 0
2989                MSC &= msc
2990                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2991            n = -1
2992            if debug: print (MSC)
2993            for i,M in enumerate(MSC):
2994                if M:
2995                    n += 1
2996                    CSI[0][i][0] = n+1
2997                    CSI[1][i][0] = 1.0
2998
2999        return list(CSI),dM,dMTP
3000       
3001    if debug: print ('super space group: '+SSGData['SSpGrp'])
3002    xyz = np.array(XYZ)%1.
3003    SGOps = copy.deepcopy(SGData['SGOps'])
3004    laue = SGData['SGLaue']
3005    siteSym = SytSym(XYZ,SGData)[0].strip()
3006    if debug: print ('siteSym: '+siteSym)
3007    SSGOps = copy.deepcopy(SSGData['SSGOps'])
3008    #expand ops to include inversions if any
3009    if SGData['SGInv'] and not SGData['SGFixed']:
3010        for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
3011            SGOps.append([-op[0],-op[1]%1.])
3012            SSGOps.append([-sop[0],-sop[1]%1.])
3013    #build set of sym ops around special position       
3014    SSop = []
3015    Sop = []
3016    Sdtau = []
3017    for iop,Op in enumerate(SGOps):         
3018        nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
3019        if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
3020            SSop.append(SSGOps[iop])
3021            Sop.append(SGOps[iop])
3022            ssopinv = nl.inv(SSGOps[iop][0])
3023            mst = ssopinv[3][:3]
3024            epsinv = ssopinv[3][3]
3025            Sdtau.append(np.sum(mst*(XYZ-SGOps[iop][1])-epsinv*SSGOps[iop][1][3]))
3026    SdIndx = np.argsort(np.array(Sdtau))     # just to do in sensible order
3027    if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
3028    #setup displacement arrays
3029    tau = np.linspace(-1,1,49,True)
3030    #make modulation arrays - one parameter at a time
3031    if Stype == 'Sfrac':
3032        CSI,dF,dFTP = DoFrac()
3033    elif Stype == 'Spos':
3034        CSI,dF,dFTP = DoXYZ()
3035        CSI[0] = orderParms(CSI[0])
3036    elif Stype == 'Sadp':
3037        CSI,dF,dFTP = DoUij()
3038        CSI[0] = orderParms(CSI[0]) 
3039    elif Stype == 'Smag':
3040        CSI,dF,dFTP = DoMag()
3041
3042    if debug:
3043        return CSI,dF,dFTP
3044    else:
3045        return CSI,[],[]
3046   
3047def MustrainNames(SGData):
3048    'Needs a doc string'
3049    laue = SGData['SGLaue']
3050    uniq = SGData['SGUniq']
3051    if laue in ['m3','m3m']:
3052        return ['S400','S220']
3053    elif laue in ['6/m','6/mmm','3m1']:
3054        return ['S400','S004','S202']
3055    elif laue in ['31m','3']:
3056        return ['S400','S004','S202','S301']
3057    elif laue in ['3R','3mR']:
3058        return ['S400','S220','S310','S211']
3059    elif laue in ['4/m','4/mmm']:
3060        return ['S400','S004','S220','S022']
3061    elif laue in ['mmm']:
3062        return ['S400','S040','S004','S220','S202','S022']
3063    elif laue in ['2/m']:
3064        SHKL = ['S400','S040','S004','S220','S202','S022']
3065        if uniq == 'a':
3066            SHKL += ['S013','S031','S211']
3067        elif uniq == 'b':
3068            SHKL += ['S301','S103','S121']
3069        elif uniq == 'c':
3070            SHKL += ['S130','S310','S112']
3071        return SHKL
3072    else:
3073        SHKL = ['S400','S040','S004','S220','S202','S022']
3074        SHKL += ['S310','S103','S031','S130','S301','S013']
3075        SHKL += ['S211','S121','S112']
3076        return SHKL
3077       
3078def HStrainVals(HSvals,SGData):
3079    laue = SGData['SGLaue']
3080    uniq = SGData['SGUniq']
3081    DIJ = np.zeros(6)
3082    if laue in ['m3','m3m']:
3083        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
3084    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3085        DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
3086    elif laue in ['3R','3mR']:
3087        DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
3088    elif laue in ['4/m','4/mmm']:
3089        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
3090    elif laue in ['mmm']:
3091        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3092    elif laue in ['2/m']:
3093        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3094        if uniq == 'a':
3095            DIJ[5] = HSvals[3]
3096        elif uniq == 'b':
3097            DIJ[4] = HSvals[3]
3098        elif uniq == 'c':
3099            DIJ[3] = HSvals[3]
3100    else:
3101        DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
3102    return DIJ
3103
3104def HStrainNames(SGData):
3105    'Needs a doc string'
3106    laue = SGData['SGLaue']
3107    uniq = SGData['SGUniq']
3108    if laue in ['m3','m3m']:
3109        return ['D11','eA']         #add cubic strain term
3110    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3111        return ['D11','D33']
3112    elif laue in ['3R','3mR']:
3113        return ['D11','D12']
3114    elif laue in ['4/m','4/mmm']:
3115        return ['D11','D33']
3116    elif laue in ['mmm']:
3117        return ['D11','D22','D33']
3118    elif laue in ['2/m']:
3119        Dij = ['D11','D22','D33']
3120        if uniq == 'a':
3121            Dij += ['D23']
3122        elif uniq == 'b':
3123            Dij += ['D13']
3124        elif uniq == 'c':
3125            Dij += ['D12']
3126        return Dij
3127    else:
3128        Dij = ['D11','D22','D33','D12','D13','D23']
3129        return Dij
3130   
3131def MustrainCoeff(HKL,SGData):
3132    'Needs a doc string'
3133    #NB: order of terms is the same as returned by MustrainNames
3134    laue = SGData['SGLaue']
3135    uniq = SGData['SGUniq']
3136    h,k,l = HKL
3137    Strm = []
3138    if laue in ['m3','m3m']:
3139        Strm.append(h**4+k**4+l**4)
3140        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
3141    elif laue in ['6/m','6/mmm','3m1']:
3142        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
3143        Strm.append(l**4)
3144        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
3145    elif laue in ['31m','3']:
3146        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
3147        Strm.append(l**4)
3148        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
3149        Strm.append(4.0*l*h**3)
3150    elif laue in ['3R','3mR']:
3151        Strm.append(h**4+k**4+l**4)
3152        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
3153        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
3154        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
3155    elif laue in ['4/m','4/mmm']:
3156        Strm.append(h**4+k**4)
3157        Strm.append(l**4)
3158        Strm.append(3.0*(h*k)**2)
3159        Strm.append(3.0*((h*l)**2+(k*l)**2))
3160    elif laue in ['mmm']:
3161        Strm.append(h**4)
3162        Strm.append(k**4)
3163        Strm.append(l**4)
3164        Strm.append(3.0*(h*k)**2)
3165        Strm.append(3.0*(h*l)**2)
3166        Strm.append(3.0*(k*l)**2)
3167    elif laue in ['2/m']:
3168        Strm.append(h**4)
3169        Strm.append(k**4)
3170        Strm.append(l**4)
3171        Strm.append(3.0*(h*k)**2)
3172        Strm.append(3.0*(h*l)**2)
3173        Strm.append(3.0*(k*l)**2)
3174        if uniq == 'a':
3175            Strm.append(2.0*k*l**3)
3176            Strm.append(2.0*l*k**3)
3177            Strm.append(4.0*k*l*h**2)
3178        elif uniq == 'b':
3179            Strm.append(2.0*l*h**3)
3180            Strm.append(2.0*h*l**3)
3181            Strm.append(4.0*h*l*k**2)
3182        elif uniq == 'c':
3183            Strm.append(2.0*h*k**3)
3184            Strm.append(2.0*k*h**3)
3185            Strm.append(4.0*h*k*l**2)
3186    else:
3187        Strm.append(h**4)
3188        Strm.append(k**4)
3189        Strm.append(l**4)
3190        Strm.append(3.0*(h*k)**2)
3191        Strm.append(3.0*(h*l)**2)
3192        Strm.append(3.0*(k*l)**2)
3193        Strm.append(2.0*k*h**3)
3194        Strm.append(2.0*h*l**3)
3195        Strm.append(2.0*l*k**3)
3196        Strm.append(2.0*h*k**3)
3197        Strm.append(2.0*l*h**3)
3198        Strm.append(2.0*k*l**3)
3199        Strm.append(4.0*k*l*h**2)
3200        Strm.append(4.0*h*l*k**2)
3201        Strm.append(4.0*k*h*l**2)
3202    return Strm
3203
3204def MuShklMean(SGData,Amat,Shkl):
3205   
3206    def genMustrain(xyz,Shkl):
3207        uvw = np.inner(Amat.T,xyz)
3208        Strm = np.array(MustrainCoeff(uvw,SGData))
3209        Sum = np.sum(np.multiply(Shkl,Strm))
3210        Sum = np.where(Sum > 0.01,Sum,0.01)
3211        Sum = np.sqrt(Sum)
3212        return Sum*xyz
3213       
3214    PHI = np.linspace(0.,360.,30,True)
3215    PSI = np.linspace(0.,180.,30,True)
3216    X = np.outer(npcosd(PHI),npsind(PSI))
3217    Y = np.outer(npsind(PHI),npsind(PSI))
3218    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3219    XYZ = np.dstack((X,Y,Z))
3220    XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
3221    return np.sqrt(np.sum(XYZ**2)/900.)
3222   
3223def Muiso2Shkl(muiso,SGData,cell):
3224    "this is to convert isotropic mustrain to generalized Shkls"
3225    import GSASIIlattice as G2lat
3226    A = G2lat.cell2AB(cell)[0]
3227   
3228    def minMus(Shkl,muiso,H,SGData,A):
3229        U = np.inner(A.T,H)
3230        S = np.array(MustrainCoeff(U,SGData))
3231        nS = S.shape[0]
3232        Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:nS,nxs]),axis=0))
3233        rad = np.sqrt(np.sum((Sum[:,nxs]*H)**2,axis=1))
3234        return (muiso-rad)**2
3235       
3236    laue = SGData['SGLaue']
3237    PHI = np.linspace(0.,360.,60,True)
3238    PSI = np.linspace(0.,180.,60,True)
3239    X = np.outer(npsind(PHI),npsind(PSI))
3240    Y = np.outer(npcosd(PHI),npsind(PSI))
3241    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3242    HKL = np.dstack((X,Y,Z))
3243    if laue in ['m3','m3m']:
3244        S0 = [1000.,1000.]
3245    elif laue in ['6/m','6/mmm']:
3246        S0 = [1000.,1000.,1000.]
3247    elif laue in ['31m','3','3m1']:
3248        S0 = [1000.,1000.,1000.,1000.]
3249    elif laue in ['3R','3mR']:
3250        S0 = [1000.,1000.,1000.,1000.]
3251    elif laue in ['4/m','4/mmm']:
3252        S0 = [1000.,1000.,1000.,1000.]
3253    elif laue in ['mmm']:
3254        S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
3255    elif laue in ['2/m']:
3256        S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
3257    else:
3258        S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000., 
3259            1000.,1000.,0.,0.,0.]
3260    S0 = np.array(S0)
3261    HKL = np.reshape(HKL,(-1,3))
3262    result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
3263    return result[0]
3264       
3265def PackRot(SGOps):
3266    IRT = []
3267    for ops in SGOps:
3268        M = ops[0]
3269        irt = 0
3270        for j in range(2,-1,-1):
3271            for k in range(2,-1,-1):
3272                irt *= 3
3273                irt += M[k][j]
3274        IRT.append(int(irt))
3275    return IRT
3276       
3277def SytSym(XYZ,SGData):
3278    '''
3279    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
3280
3281    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
3282    :param SGData: from SpcGroup
3283    :Returns: a four element tuple:
3284
3285     * The 1st element is a code for the site symmetry (see GetKNsym)
3286     * The 2nd element is the site multiplicity
3287     * Ndup number of overlapping operators
3288     * dupDir Dict - dictionary of overlapping operators
3289
3290    '''
3291    Mult = 1
3292    Isym = 0
3293    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
3294        Isym = 1073741824
3295    Jdup = 0
3296    Ndup = 0
3297    dupDir = {}
3298    inv = SGData['SGInv']+1
3299    icen = SGData['SGCen']
3300    Ncen = len(icen)
3301    if SGData['SGFixed']:       #already in list of operators
3302        inv = 1
3303    if SGData['SGGray'] and Ncen > 1: Ncen //= 2
3304    Xeqv = list(GenAtom(XYZ,SGData,True))
3305#    for xeqv in Xeqv:   print(xeqv)
3306    IRT = PackRot(SGData['SGOps'])
3307    L = -1
3308    for ic,cen in enumerate(icen[:Ncen]):
3309        for invers in range(int(inv)):
3310            for io,ops in enumerate(SGData['SGOps']):
3311                irtx = (1-2*invers)*IRT[io]
3312                L += 1
3313                if not Xeqv[L][1]:
3314                    Ndup = io
3315                    Jdup += 1
3316                    jx = GetOprPtrNumber(str(irtx))   #[KN table no,op name,KNsym ptr]
3317                    if jx < 39:
3318                        px = GetOprName(str(irtx))
3319                        if Xeqv[L][-1] < 0:
3320                            if '(' in px:
3321                                px = px.split('(')
3322                                px[0] += "'"
3323                                px = '('.join(px)
3324                            else:   
3325                                px += "'"
3326                        dupDir[px] = L
3327                        Isym += 2**(jx-1)
3328    if Isym == 1073741824: Isym = 0
3329    Mult = len(SGData['SGOps'])*Ncen*inv//Jdup
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 == '1':       #genersl 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,Mom):
3570    'Needs a doc string'
3571    SGOps = SGData['SGOps']
3572    Ax = A.split('+')
3573    Ax[0] = int(Ax[0])
3574    iAx = abs(Ax[0])
3575    nA = iAx%100-1
3576    if SGData['SGInv'] and not SGData['SGFixed']:
3577        nC = 2*len(SGOps)*(iAx//100)
3578    else:
3579        nC = len(SGOps)*(iAx//100)
3580    NA = nA
3581    if Ax[0] < 0:
3582        NA += len(SGOps)
3583    M,T = SGOps[nA]
3584    if SGData['SGGray']:
3585        newMom = -np.inner(Mom,M).T*nl.det(M)*SGData['SpnFlp'][NA+nC]
3586    else:
3587        newMom = np.inner(Mom,M).T*nl.det(M)*SGData['SpnFlp'][NA+nC]
3588#        print(len(SGOps),Ax[0],iAx,nC,nA,NA,MT2text([M,T]).replace(' ',''),SGData['SpnFlp'][NA],Mom,newMom)
3589#    print(Mom,newMom,MT2text([M,T]),)
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    """
3603    SGOps = SGData['SGOps']