source: trunk/GSASIIspc.py @ 4810

Last change on this file since 4810 was 4810, checked in by vondreele, 10 months ago

fix depreciation warning in G2indx.findMV
fix numpy round bug in G2lattice.MaxIndex?
fix bug in G2pwdGUI.OnFindMV
fix bug in G2spc.SSpcGroup

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 174.9 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: 2021-02-12 15:51:39 +0000 (Fri, 12 Feb 2021) $
12# $Author: vondreele $
13# $Revision: 4810 $
14# $URL: trunk/GSASIIspc.py $
15# $Id: GSASIIspc.py 4810 2021-02-12 15:51:39Z vondreele $
16########### SVN repository information ###################
17from __future__ import division, print_function
18import numpy as np
19import numpy.linalg as nl
20import scipy.optimize as so
21import sys
22import copy
23import os.path as ospath
24
25import GSASIIpath
26GSASIIpath.SetVersionNumber("$Revision: 4810 $")
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            if mod:
1405                return [-SSGKl[i] if mod[i] in ['a','b','g'] else SSGKl[i] for i in range(3)]
1406            else:
1407                return [SSGKl[i] for i in range(3)]
1408       
1409    def extendSSGOps(SSGOps):
1410        for OpA in SSGOps:
1411            OpAtxt = SSMT2text(OpA)
1412            if 't' not in OpAtxt:
1413                continue
1414            for OpB in SSGOps:
1415                OpBtxt = SSMT2text(OpB)
1416                if 't' not in OpBtxt:
1417                    continue
1418                OpC = list(SGProd(OpB,OpA))
1419                OpC[1] %= 1.
1420                OpCtxt = SSMT2text(OpC)
1421#                print OpAtxt.replace(' ','')+' * '+OpBtxt.replace(' ','')+' = '+OpCtxt.replace(' ','')
1422                for k,OpD in enumerate(SSGOps):
1423                    OpDtxt = SSMT2text(OpD)
1424                    OpDtxt2 = ''
1425                    if SGData['SGGray']:                       
1426                        OpDtxt2 = SSMT2text([OpD[0],OpD[1]+np.array([0.,0.,0.,.5])])
1427#                    print '    ('+OpCtxt.replace(' ','')+' = ? '+OpDtxt.replace(' ','')+')'
1428                    if OpCtxt == OpDtxt:
1429                        continue
1430                    elif OpCtxt == OpDtxt2:
1431                        continue
1432                    elif OpCtxt.split(',')[:3] == OpDtxt.split(',')[:3]:
1433                        if 't' not in OpDtxt:
1434                            SSGOps[k] = OpC
1435#                            print k,'   new:',OpCtxt.replace(' ','')
1436                            break
1437                        else:
1438                            OpCtxt = OpCtxt.replace(' ','')
1439                            OpDtxt = OpDtxt.replace(' ','')
1440                            Txt = OpCtxt+' conflicts with '+OpDtxt
1441#                            print (Txt)
1442                            return False,Txt
1443        return True,SSGOps
1444       
1445    def findMod(modSym):
1446        for a in ['a','b','g']:
1447            if a in modSym:
1448                return a
1449               
1450    def genSSGOps():
1451        SSGOps = SSGData['SSGOps'][:]
1452        iFrac = {}
1453        for i,frac in enumerate(SSGData['modSymb']):
1454            if frac in ['1/2','1/3','1/4','1/6','1']:
1455                iFrac[i] = frac+'.'
1456#        print SGData['SpGrp']+SSymbol
1457#        print 'SSGKl',SSGKl,'genQ',genQ,'iFrac',iFrac,'modSymb',SSGData['modSymb']
1458# set identity & 1,-1; triclinic
1459        SSGOps[0][0][3,3] = 1.
1460## expand if centrosymmetric
1461#        if SGData['SGInv']:
1462#            SSGOps += [[-1*M,V] for M,V in SSGOps[:]]
1463# monoclinic - all done & all checked
1464        if SGData['SGPtGrp'] in ['2','m']:  #OK
1465            SSGOps[1][0][3,3] = SSGKl[0]
1466            SSGOps[1][1][3] = genQ[0]
1467            for i in iFrac:
1468                SSGOps[1][0][3,i] = -SSGKl[0]
1469        elif SGData['SGPtGrp'] == '2/m':    #OK
1470            SSGOps[1][0][3,3] = SSGKl[1]
1471            if 's' in gensym:
1472                SSGOps[1][1][3] = 0.5
1473            for i in iFrac:
1474                SSGOps[1][0][3,i] = SSGKl[0]
1475           
1476# orthorhombic - all OK not fully checked
1477        elif SGData['SGPtGrp'] in ['222','mm2','m2m','2mm']:    #OK
1478            if SGData['SGPtGrp'] == '222':
1479                OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[1,2],2:[1,3]},'b':{2:[3,2],0:[1,2]}} #OK
1480            elif SGData['SGPtGrp'] == 'mm2':
1481                OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]}} #OK
1482            elif SGData['SGPtGrp'] == 'm2m':
1483                OrOps = {'b':{0:[1,2],2:[3,2]},'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]}} #OK
1484            elif SGData['SGPtGrp'] == '2mm':
1485                OrOps = {'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]},'g':{0:[1,3],1:[2,3]}} #OK
1486            a = findMod(SSGData['modSymb'])
1487            OrFrac = OrOps[a]
1488            for j in iFrac:
1489                for i in OrFrac[j]:
1490                    SSGOps[i][0][3,j] = -2.*eval(iFrac[j])*SSGKl[i-1]
1491            for i in [0,1,2]:
1492                SSGOps[i+1][0][3,3] = SSGKl[i]
1493                SSGOps[i+1][1][3] = genQ[i]
1494                E,SSGOps = extendSSGOps(SSGOps)
1495                if not E:
1496                    return E,SSGOps
1497        elif SGData['SGPtGrp'] == 'mmm':    #OK
1498            OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]}} 
1499            a = findMod(SSGData['modSymb'])
1500            if a == 'g':
1501                SSkl = [1,1,1]
1502            elif a == 'a':
1503                SSkl = [-1,1,-1]
1504            else:
1505                SSkl = [1,-1,-1]
1506            OrFrac = OrOps[a]
1507            for j in iFrac:
1508                for i in OrFrac[j]:
1509                    SSGOps[i][0][3,j] = -2.*eval(iFrac[j])*SSkl[i-1]
1510            for i in [0,1,2]:
1511                SSGOps[i+1][0][3,3] = SSkl[i]
1512                SSGOps[i+1][1][3] = genQ[i]
1513                E,SSGOps = extendSSGOps(SSGOps)
1514                if not E:
1515                    return E,SSGOps               
1516# tetragonal - all done & checked
1517        elif SGData['SGPtGrp'] == '4':  #OK
1518            SSGOps[1][0][3,3] = SSGKl[0]
1519            SSGOps[1][1][3] = genQ[0]
1520            if '1/2' in SSGData['modSymb']:
1521                SSGOps[1][0][3,1] = -1
1522        elif SGData['SGPtGrp'] == '-4': #OK
1523            SSGOps[1][0][3,3] = SSGKl[0]
1524            if '1/2' in SSGData['modSymb']:
1525                SSGOps[1][0][3,1] = 1
1526        elif SGData['SGPtGrp'] in ['4/m',]: #OK
1527            if '1/2' in SSGData['modSymb']:
1528                SSGOps[1][0][3,1] = -SSGKl[0]
1529            for i,j in enumerate([1,3]):
1530                SSGOps[j][0][3,3] = 1
1531                if genQ[i]:
1532                    SSGOps[j][1][3] = genQ[i]
1533                E,SSGOps = extendSSGOps(SSGOps)
1534                if not E:
1535                    return E,SSGOps
1536        elif SGData['SGPtGrp'] in ['422','4mm','-42m','-4m2',]: #OK
1537            iGens = [1,4,5]
1538            if SGData['SGPtGrp'] in ['4mm','-4m2',]:
1539                iGens = [1,6,7]
1540            for i,j in enumerate(iGens):
1541                if '1/2' in SSGData['modSymb'] and i < 2:
1542                    SSGOps[j][0][3,1] = SSGKl[i]
1543                SSGOps[j][0][3,3] = SSGKl[i]
1544                if genQ[i]:
1545                    if 's' in gensym and j == 6:
1546                        SSGOps[j][1][3] = -genQ[i]
1547                    else:
1548                        SSGOps[j][1][3] = genQ[i]
1549                E,SSGOps = extendSSGOps(SSGOps)
1550                if not E:
1551                    return E,SSGOps
1552        elif SGData['SGPtGrp'] in ['4/mmm',]:#OK
1553            if '1/2' in SSGData['modSymb']:
1554                SSGOps[1][0][3,1] = -SSGKl[0]
1555                SSGOps[6][0][3,1] = SSGKl[1]
1556                if modsym:
1557                   SSGOps[1][1][3]  = -genQ[3]
1558            for i,j in enumerate([1,2,6,7]):
1559                SSGOps[j][0][3,3] = 1
1560                SSGOps[j][1][3] = genQ[i]
1561                E,Result = extendSSGOps(SSGOps)
1562                if not E:
1563                    return E,Result
1564                else:
1565                    SSGOps = Result
1566               
1567# trigonal - all done & checked
1568        elif SGData['SGPtGrp'] == '3':  #OK
1569            SSGOps[1][0][3,3] = SSGKl[0]
1570            if '1/3' in SSGData['modSymb']:
1571                SSGOps[1][0][3,1] = -1
1572            SSGOps[1][1][3] = genQ[0]
1573        elif SGData['SGPtGrp'] == '-3': #OK
1574            SSGOps[1][0][3,3] = -SSGKl[0]
1575            if '1/3' in SSGData['modSymb']:
1576                SSGOps[1][0][3,1] = -1
1577            SSGOps[1][1][3] = genQ[0]
1578        elif SGData['SGPtGrp'] in ['312','3m','-3m','-3m1','3m1']:   #OK
1579            if '1/3' in SSGData['modSymb']:
1580                SSGOps[1][0][3,1] = -1
1581            for i,j in enumerate([1,5]):
1582                if SGData['SGPtGrp'] in ['3m','-3m']:
1583                    SSGOps[j][0][3,3] = 1
1584                else:                   
1585                    SSGOps[j][0][3,3] = SSGKl[i+1]
1586                if genQ[i]:
1587                    SSGOps[j][1][3] = genQ[i]
1588        elif SGData['SGPtGrp'] in ['321','32']:   #OK
1589            for i,j in enumerate([1,4]):
1590                SSGOps[j][0][3,3] = SSGKl[i]
1591                if genQ[i]:
1592                    SSGOps[j][1][3] = genQ[i]
1593        elif SGData['SGPtGrp'] in ['31m','-31m']:   #OK
1594            ids = [1,3]
1595            if SGData['SGPtGrp'] == '-31m':
1596                ids = [1,3]
1597            if '1/3' in SSGData['modSymb']:
1598                SSGOps[ids[0]][0][3,1] = -SSGKl[0]
1599            for i,j in enumerate(ids):
1600                SSGOps[j][0][3,3] = 1
1601                if genQ[i+1]:
1602                    SSGOps[j][1][3] = genQ[i+1]
1603                     
1604# hexagonal all done & checked
1605        elif SGData['SGPtGrp'] == '6':  #OK
1606            SSGOps[1][0][3,3] = SSGKl[0]
1607            SSGOps[1][1][3] = genQ[0]
1608        elif SGData['SGPtGrp'] == '-6': #OK
1609            SSGOps[1][0][3,3] = SSGKl[0]
1610        elif SGData['SGPtGrp'] in ['6/m',]: #OK
1611            SSGOps[1][0][3,3] = -SSGKl[1]
1612            SSGOps[1][1][3] = genQ[0]
1613            SSGOps[2][1][3] = genQ[1]
1614        elif SGData['SGPtGrp'] in ['622',]: #OK
1615            for i,j in enumerate([1,9,8]):
1616                SSGOps[j][0][3,3] = SSGKl[i]
1617                if genQ[i]:
1618                    SSGOps[j][1][3] = -genQ[i]
1619                E,SSGOps = extendSSGOps(SSGOps)
1620           
1621        elif SGData['SGPtGrp'] in ['6mm','-62m','-6m2',]: #OK
1622            for i,j in enumerate([1,6,7]):
1623                SSGOps[j][0][3,3] = SSGKl[i]
1624                if genQ[i]:
1625                    SSGOps[j][1][3] = genQ[i]
1626                E,SSGOps = extendSSGOps(SSGOps)
1627        elif SGData['SGPtGrp'] in ['6/mmm',]: # OK
1628            for i,j in enumerate([1,2,10,11]):
1629                SSGOps[j][0][3,3] = 1
1630                if genQ[i]:
1631                    SSGOps[j][1][3] = genQ[i]
1632                E,SSGOps = extendSSGOps(SSGOps)
1633        elif SGData['SGPtGrp'] in ['1','-1']: #triclinic - done
1634            return True,SSGOps
1635        E,SSGOps = extendSSGOps(SSGOps)
1636        return E,SSGOps
1637       
1638    def specialGen(gensym,modsym):
1639        sym = ''.join(gensym)
1640        if SGData['SGPtGrp'] in ['2/m',] and 'n' in SGData['SpGrp']:
1641            if 's' in sym:
1642                gensym = 'ss'
1643        if SGData['SGPtGrp'] in ['-62m',] and sym == '00s':
1644            gensym = '0ss'
1645        elif SGData['SGPtGrp'] in ['222',]:
1646            if sym == '00s':
1647                gensym = '0ss'
1648            elif sym == '0s0':
1649                gensym = 'ss0'
1650            elif sym == 's00':
1651                gensym = 's0s'
1652        elif SGData['SGPtGrp'] in ['mmm',]:
1653            if 'g' in modsym:
1654                if sym == 's00':
1655                    gensym = 's0s'
1656                elif sym == '0s0':
1657                    gensym = '0ss'
1658            elif 'a' in modsym:
1659                if sym == '0s0':
1660                    gensym = 'ss0'
1661                elif sym == '00s':
1662                    gensym = 's0s'
1663            elif 'b' in modsym:
1664                if sym == '00s':
1665                    gensym = '0ss'
1666                elif sym == 's00':
1667                    gensym = 'ss0'
1668        return gensym
1669                           
1670    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.}
1671    if SGData['SGLaue'] in ['m3','m3m']:
1672        return '(3+1) superlattices not defined for cubic space groups',None
1673    elif SGData['SGLaue'] in ['3R','3mR']:
1674        return '(3+1) superlattices not defined for rhombohedral settings - use hexagonal setting',None
1675    try:
1676        modsym,gensym = splitSSsym(SSymbol)
1677    except ValueError:
1678        return 'Error in superspace symbol '+SSymbol,None
1679    modQ = [Fracs[mod] for mod in modsym]
1680    SSGKl = SGData['SSGKl'][:]
1681    if SGData['SGLaue'] in ['2/m','mmm']:
1682        SSGKl = fixMonoOrtho()
1683    Ngen = len(gensym)
1684    if SGData.get('SGGray',False):
1685        Ngen -= 1
1686    if len(gensym) and Ngen != len(SSGKl):
1687        return 'Wrong number of items in generator symbol '+''.join(gensym),None
1688    gensym = specialGen(gensym[:Ngen],modsym)
1689    genQ = [Fracs[mod] for mod in gensym[:Ngen]]
1690    if not genQ:
1691        genQ = [0,0,0,0]
1692    SSgSpc = SGData['SpGrp']+SSymbol
1693    if SGData['SGGray']:
1694        SSgSpc = SSgSpc.replace('('," 1'(")
1695    SSGData = {'SSpGrp':SSgSpc,'modQ':modQ,'modSymb':modsym,'SSGKl':SSGKl}
1696    SSCen = np.zeros((len(SGData['SGCen']),4))
1697    for icen,cen in enumerate(SGData['SGCen']):
1698        SSCen[icen,0:3] = cen
1699    if 'BNSlattsym' in SGData and '_' in SGData['BNSlattsym'][0]:
1700        Ncen = len(SGData['SGCen'])
1701        for icen in range(Ncen//2,Ncen):
1702            SSCen[icen,3] = 0.5
1703    SSGData['SSGCen'] = SSCen%1.
1704    SSGData['SSGOps'] = []
1705    for iop,op in enumerate(SGData['SGOps']):
1706        T = np.zeros(4)
1707        ssop = np.zeros((4,4))
1708        ssop[:3,:3] = op[0]
1709        T[:3] = op[1]
1710        SSGData['SSGOps'].append([ssop,T])
1711    E,Result = genSSGOps()
1712    if E:
1713        SSGData['SSGOps'] = Result
1714        if DEBUG:
1715            print ('Super spacegroup operators for '+SSGData['SSpGrp'])
1716            for Op in Result:
1717                print (SSMT2text(Op).replace(' ',''))
1718            if SGData['SGInv']:                                 
1719                for Op in Result:
1720                    Op = [-Op[0],-Op[1]%1.]
1721                    print (SSMT2text(Op).replace(' ',''))                                 
1722        return None,SSGData
1723    else:
1724        return Result+'\nOperator conflict - incorrect superspace symbol',None
1725   
1726def SSChoice(SGData):
1727    '''
1728    Gets the unique set of possible super space groups for a given space group
1729    '''
1730    ptgpSS = {'1':['(abg)',],'-1':['(abg)',],
1731                   
1732        '2':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1733        'm':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1734        '2/m':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1735       
1736        '222':['(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        'mm2':['(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        'm2m':['(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        '2mm':['(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        'mmm':['(00g)','(1/20g)','(01/2g)','(1/21/2g)','(10g)','(01g)',
1749               '(a00)','(a1/20)','(a01/2)','(a1/21/2)','(a10)','(a01)',
1750               '(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)','(1b0)','(0b1)',],
1751               
1752        '4':['(00g)','(1/21/2g)'],'4mm':['(00g)','(1/21/2g)'],
1753        '4/m':['(00g)','(1/21/2g)'],
1754        '422':['(00g)','(1/21/2g)'],'-4m2':['(00g)','(1/21/2g)'],'-42m':['(00g)','(1/21/2g)'],
1755        '4/mmm':['(00g)','(1/21/2g)'],
1756       
1757        '3':['(00g)','(1/31/3g)'],'-3':['(00g)','(1/31/3g)'],
1758        '32':['(00g)'],'3m':['(00g)'],'-3m':['(00g)'],
1759        '321':['(00g)'],'3m1':['(00g)'],'-3m1':['(00g)'],
1760        '312':['(00g)','(1/31/3g)'],'31m':['(00g)','(1/31/3g)'],'-31m':['(00g)','(1/31/3g)'],
1761       
1762        '6':['(00g)',],'6/m':['(00g)',],'-62m':['(00g)',],'-6m2':['(00g)',],
1763        '622':['(00g)',],'6/mmm':['(00g)',],'6mm':['(00g)',],
1764       
1765        '23':['',],'m3':['',],'432':['',],'-43m':['',],'m3m':['',]}
1766           
1767    ptgpTS = {'1':['0',],'-1':['0',],
1768             
1769        '2':['0','s'],'m':['0','s'],
1770        '2/m':['00','0s','ss','s0'],
1771       
1772        '222':['000','s00','0s0','00s',],
1773        'mm2':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1774        'm2m':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1775        '2mm':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1776        'mmm':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1777       
1778        '4':['0','q','s'],'4mm':['000','q00','s00','s0s','ss0','0ss','qq0','qqs'],
1779        '4/m':['00','s0'],'-4m2':['000','0s0','0q0'],'-42m':['000','00s'],
1780        '422':['000','q00','s00','s0s','ss0','0ss','qq0','qqs','0q0'],
1781        '4/mmm':['0000','s0s0','00ss','s00s','ss00','0ss0','0s0s'],
1782       
1783        '3':['0','t'],'-3':['0','t'],
1784        '32':['00','t0'],'3m':['00','0s'],'-3m':['00','0s'],
1785        '321':['000','t00'],'3m1':['000','0s0'],'-3m1':['000','0s0'],
1786        '312':['000','t00'],'31m':['000','00s'],'-31m':['000','00s'],
1787       
1788        '6':['0','h','t','s'],
1789        '6/m':['00','s0'],'-62m':['000','00s'],'-6m2':['000','0s0'],
1790        '622':['000','h00','t00','s00',],'6mm':['000','ss0','s0s','0ss',],
1791        '6/mmm':['0000','s0s0','00ss','s00s','ss00','0ss0','0s0s'],
1792       
1793        '23':['',],'m3':['',],'432':['',],'-43m':['',],'m3m':['',]}
1794   
1795    ptgp = SGData['SGPtGrp']
1796    SSChoice = []
1797    for ax in ptgpSS[ptgp]:
1798        for sx in ptgpTS[ptgp]:
1799            SSChoice.append(ax+sx)
1800            if SGData['SGGray']: SSChoice[-1] += 's'
1801    ssChoice = []
1802    ssHash = []
1803    for item in SSChoice:
1804        E,SSG = SSpcGroup(SGData,item)
1805        if SSG:
1806            sshash = hash(str(SSGPrint(SGData,SSG)[1]))
1807            if sshash not in ssHash:
1808                ssHash.append(sshash)
1809                ssChoice.append(item)
1810    return ssChoice
1811           
1812def splitSSsym(SSymbol):
1813    '''
1814    Splits supersymmetry symbol into two lists of strings
1815    '''
1816    mssym = SSymbol.replace(' ','').split(')')
1817    if len(mssym) > 1:
1818        modsym,gensym = mssym
1819    else:
1820        modsym = mssym[0]
1821        gensym = ''
1822    modsym = modsym.replace(',','')
1823    if "1'" in modsym:
1824        gensym = gensym[:-1]
1825    modsym = modsym.replace("1'",'')
1826    if gensym in ['0','00','000','0000']:       #get rid of extraneous symbols
1827        gensym = ''
1828    nfrac = modsym.count('/')
1829    modsym = modsym.lstrip('(')
1830    if nfrac == 0:
1831        modsym = list(modsym)
1832    elif nfrac == 1:
1833        pos = modsym.find('/')
1834        if pos == 1:
1835            modsym = [modsym[:3],modsym[3],modsym[4]]
1836        elif pos == 2:
1837            modsym = [modsym[0],modsym[1:4],modsym[4]]
1838        else:
1839            modsym = [modsym[0],modsym[1],modsym[2:]]
1840    else:
1841        lpos = modsym.find('/')
1842        rpos = modsym.rfind('/')
1843        if lpos == 1 and rpos == 4:
1844            modsym = [modsym[:3],modsym[3:6],modsym[6]]
1845        elif lpos == 1 and rpos == 5:
1846            modsym = [modsym[:3],modsym[3],modsym[4:]]
1847        else:
1848            modsym = [modsym[0],modsym[1:4],modsym[4:]]
1849    gensym = list(gensym)
1850    return modsym,gensym
1851       
1852def SSGPrint(SGData,SSGData,AddInv=False):
1853    '''
1854    Print the output of SSpcGroup in a nicely formatted way. Used in SSpaceGroup
1855
1856    :param SGData: space group data structure as defined in SpcGroup above.
1857    :param SSGData: from :func:`SSpcGroup`
1858    :returns:
1859        SSGText - list of strings with the superspace group details
1860        SGTable - list of strings for each of the operations
1861    '''
1862    nCen = len(SSGData['SSGCen'])
1863    Mult = nCen*len(SSGData['SSGOps'])*(int(SGData['SGInv'])+1)
1864    if SGData.get('SGFixed',False):
1865        Mult = len(SSGData['SSGCen'])*len(SSGData['SSGOps'])
1866    SSsymb = SSGData['SSpGrp']
1867    if 'BNSlattsym' in SGData and '_' in SGData['BNSlattsym'][0]:
1868        SSsymb = SGData['BNSlattsym'][0]+SSsymb[1:]
1869    SSGCen = list(SSGData['SSGCen'])
1870    if SGData.get('SGGray',False):
1871        if SGData.get('SGFixed',False): 
1872            Mult //= 2
1873        else:
1874            SSGCen += list(SSGData['SSGCen']+[0,0,0,0.5])
1875            SSGCen =  np.array(SSGCen)%1.
1876    else:
1877        if "1'" in SSsymb:  #leftover in nonmag phase in mcif file
1878            nCen //= 2
1879            Mult //= 2
1880            SSsymb = SSsymb.replace("1'",'')[:-1]
1881    SSGText = []
1882    SSGText.append(' Superspace Group: '+SSsymb)
1883    CentStr = 'centrosymmetric'
1884    if not SGData['SGInv']:
1885        CentStr = 'non'+CentStr
1886    if SGData['SGLatt'] in 'ABCIFR':
1887        SSGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
1888    else:
1889        SSGText.append(' The superlattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())
1890    SSGText.append(' The Laue symmetry is '+SGData['SGLaue'])
1891    SGptGp = SGData['SGPtGrp']
1892    if SGData['SGGray']:
1893        SGptGp += "1'"
1894    SSGText.append(' The superlattice point group is '+SGptGp+', '+''.join([str(i) for i in SSGData['SSGKl']]))
1895    SSGText.append(' The number of superspace group generators is '+str(len(SGData['SSGKl'])))
1896    SSGText.append(' Multiplicity of a general site is '+str(Mult))
1897    if SGData['SGUniq'] in ['a','b','c']:
1898        SSGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
1899    if SGData['SGInv']:
1900        SSGText.append(' The inversion center is located at 0,0,0')
1901    if SGData['SGPolax']:
1902        SSGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
1903    SSGText.append(' ')
1904    if len(SSGCen) > 1:
1905        SSGText.append(' The equivalent positions are:')
1906        SSGText.append(' ('+SSLatt2text(SSGCen)+')+\n')
1907    else:
1908        SSGText.append(' The equivalent positions are:\n')
1909    SSGTable = []
1910    for i,Opr in enumerate(SSGData['SSGOps']):
1911        SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr)))
1912    if AddInv and SGData['SGInv']:
1913        for i,Opr in enumerate(SSGData['SSGOps']):
1914            IOpr = [-Opr[0],-Opr[1]]
1915            SSGTable.append('(%2d) %s'%(i+1+len(SSGData['SSGOps']),SSMT2text(IOpr)))       
1916    return SSGText,SSGTable
1917   
1918def SSGModCheck(Vec,modSymb,newMod=True):
1919    ''' Checks modulation vector compatibility with supersymmetry space group symbol.
1920    if newMod: Superspace group symbol takes precidence & the vector will be modified accordingly
1921    '''
1922    Fracs = {'1/2':0.5,'1/3':1./3,'1':1.0,'0':0.,'a':0.,'b':0.,'g':0.}
1923    modQ = [Fracs[mod] for mod in modSymb]
1924    if newMod:
1925        newVec = Vec
1926        if not np.any(Vec):
1927            newVec = [0.1 if (vec == 0.0 and mod in ['a','b','g']) else vec for [vec,mod] in zip(Vec,modSymb)]
1928        return [Q if mod not in ['a','b','g'] and vec != Q else vec for [vec,mod,Q] in zip(newVec,modSymb,modQ)],  \
1929            [True if mod in ['a','b','g'] else False for mod in modSymb]
1930    else:
1931        return Vec,[True if mod in ['a','b','g'] else False for mod in modSymb]
1932
1933def SSMT2text(Opr):
1934    "From superspace group matrix/translation operator returns text version"
1935    XYZS = ('x','y','z','t')    #Stokes, Campbell & van Smaalen notation
1936    TRA = ('   ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
1937    Fld = ''
1938    M,T = Opr
1939    for j in range(4):
1940        IJ = ''
1941        for k in range(4):
1942            txt = str(int(round(M[j][k])))
1943            txt = txt.replace('1',XYZS[k]).replace('0','')
1944            if '2' in txt:
1945                txt += XYZS[k]
1946            if IJ and M[j][k] > 0:
1947                IJ += '+'+txt
1948            else:
1949                IJ += txt
1950        IK = int(round(T[j]*12))%12
1951        if IK:
1952            if not IJ:
1953                break
1954            if IJ[0] == '-':
1955                Fld += (TRA[IK]+IJ).rjust(8)
1956            else:
1957                Fld += (TRA[IK]+'+'+IJ).rjust(8)
1958        else:
1959            Fld += IJ.rjust(8)
1960        if j != 3: Fld += ', '
1961    return Fld
1962   
1963def SSLatt2text(SSGCen):
1964    "Lattice centering vectors to text"
1965    lattTxt = ''
1966    lattDir = {4:'1/3',6:'1/2',8:'2/3',0:'0'}
1967    for vec in SSGCen:
1968        lattTxt += ' '
1969        for item in vec:
1970            lattTxt += '%s,'%(lattDir[int(item*12)])
1971        lattTxt = lattTxt.rstrip(',')
1972        lattTxt += ';'
1973    lattTxt = lattTxt.rstrip(';').lstrip(' ')
1974    return lattTxt
1975       
1976def SSpaceGroup(SGSymbol,SSymbol):
1977    '''
1978    Print the output of SSpcGroup in a nicely formatted way.
1979
1980    :param SGSymbol: space group symbol with spaces between axial fields.
1981    :param SSymbol: superspace group symbol extension (string).
1982    :returns: nothing
1983    '''
1984
1985    E,A = SpcGroup(SGSymbol)
1986    if E > 0:
1987        print (SGErrors(E))
1988        return
1989    E,B = SSpcGroup(A,SSymbol)   
1990    if E > 0:
1991        print (E)
1992        return
1993    for l in SSGPrint(B):
1994        print (l)
1995       
1996def SGProd(OpA,OpB):
1997    '''
1998    Form space group operator product. OpA & OpB are [M,V] pairs;
1999        both must be of same dimension (3 or 4). Returns [M,V] pair
2000    '''
2001    A,U = OpA
2002    B,V = OpB
2003    M = np.inner(B,A.T)
2004    W = np.inner(B,U)+V
2005    return M,W
2006       
2007def GetLittleGrpOps(SGData,vec):
2008    ''' Find rotation part of operators that leave vec unchanged
2009   
2010    :param SGData: space group data structure as defined in SpcGroup above.
2011    :param vec: a numpy array of fractional vector coordinates
2012    :returns: Little - list of operators [M,T] that form the little gropu
2013    '''
2014    Little = []
2015    Ops = SGData['SGOps'][:]
2016    if SGData['SGInv']:
2017        Ops += [[-M,-T] for [M,T] in Ops]
2018    for [M,T] in Ops:
2019        tvec = np.inner(M,vec)%1.
2020        if np.allclose(tvec,vec%1.):
2021            Little.append([M,T])
2022    return Little
2023       
2024def MoveToUnitCell(xyz):
2025    '''
2026    Translates a set of coordinates so that all values are >=0 and < 1
2027
2028    :param xyz: a list or numpy array of fractional coordinates
2029    :returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
2030    '''
2031    XYZ = np.array(xyz)%1.
2032    cell = np.asarray(np.rint(XYZ-xyz),dtype=np.int32)
2033    return XYZ,cell
2034       
2035def Opposite(XYZ,toler=0.0002):
2036    '''
2037    Gives opposite corner, edge or face of unit cell for position within tolerance.
2038        Result may be just outside the cell within tolerance
2039
2040    :param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
2041    :param toler: unit cell fraction tolerance making opposite
2042    :returns:
2043        XYZ: dict of opposite positions; key=unit cell & always contains XYZ
2044    '''
2045    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]]
2046    TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
2047    perm = TB*perm3
2048    cperm = ['%d,%d,%d'%(i,j,k) for i,j,k in perm]
2049    D = dict(zip(cperm,perm))
2050    new = {}
2051    for key in D:
2052        new[key] = np.array(D[key])+np.array(XYZ)
2053    return new
2054       
2055def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
2056    '''
2057    Generates the equivalent positions for a specified coordinate and space group
2058
2059    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
2060    :param SGData: from :func:`SpcGroup`
2061    :param All: True return all equivalent positions including duplicates;
2062      False return only unique positions
2063    :param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
2064    :param Move: True move generated atom positions to be inside cell
2065      False do not move atoms       
2066    :return: [[XYZEquiv],Idup,[UijEquiv],spnflp]
2067
2068      *  [XYZEquiv] is list of equivalent positions (XYZ is first entry)
2069      *  Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
2070      * is centering operator number (1-4) and - is for inversion
2071        Cell = unit cell translations needed to put new positions inside cell
2072        [UijEquiv] - equivalent Uij; absent if no Uij given
2073      * +1/-1 for spin inversion of operator - empty if not magnetic
2074       
2075    '''
2076    XYZEquiv = []
2077    UijEquiv = []
2078    Idup = []
2079    Cell = []
2080    inv = int(SGData['SGInv']+1)
2081    icen = SGData['SGCen']
2082    if SGData.get('SGFixed',False):
2083        inv = 1
2084    SpnFlp = SGData.get('SpnFlp',[])
2085    spnflp = []
2086    X = np.array(XYZ)
2087    mj = 0
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)
2102                if Move:
2103                    newX,cell = MoveToUnitCell(XT)
2104                else:
2105                    newX = XT
2106                if All:
2107                    if np.allclose(newX,X,atol=0.0002):     #do we want %1. here?
2108                        idup = False
2109                else:
2110                    if True in [np.allclose(newX%1.,oldX%1.,atol=0.0002) for oldX in XYZEquiv]:
2111                        idup = False
2112                if All or idup:
2113                    XYZEquiv.append(newX)
2114                    Idup.append(idup)
2115                    Cell.append(cell)
2116                    if len(Uij):
2117                        UijEquiv.append(newUij)
2118                    if len(SpnFlp):
2119                        spnflp.append(SpnFlp[mj])
2120                    else:
2121                        spnflp.append(1)
2122                mj += 1
2123    if len(Uij):
2124        return zip(XYZEquiv,UijEquiv,Idup,Cell,spnflp)
2125    else:
2126        return zip(XYZEquiv,Idup,Cell,spnflp)
2127       
2128def GenHKL(HKL,SGData):
2129    ''' Generates all equivlent reflections including Friedel pairs
2130    :param HKL:  [h,k,l] must be integral values
2131    :param SGData: space group data obtained from SpcGroup
2132    :returns: array Uniq: equivalent reflections
2133    '''
2134   
2135    Ops = SGData['SGOps']
2136    OpM = np.array([op[0] for op in Ops])
2137    Uniq = np.inner(OpM,HKL)
2138    Uniq = list(Uniq)+list(-1*Uniq)
2139    return np.array(Uniq)
2140
2141def GenHKLf(HKL,SGData):
2142    '''
2143    Uses old GSAS Fortran routine genhkl.for
2144
2145    :param HKL:  [h,k,l] must be integral values for genhkl.for to work
2146    :param SGData: space group data obtained from SpcGroup
2147    :returns: iabsnt,mulp,Uniq,phi
2148
2149     *   iabsnt = True if reflection is forbidden by symmetry
2150     *   mulp = reflection multiplicity including Friedel pairs
2151     *   Uniq = numpy array of equivalent hkl in descending order of h,k,l
2152     *   phi = phase offset for each equivalent h,k,l
2153
2154    '''
2155    hklf = list(HKL)+[0,]       #could be numpy array!
2156    Ops = SGData['SGOps']
2157    OpM = np.array([op[0] for op in Ops],order='F')
2158    OpT = np.array([op[1] for op in Ops])
2159    Cen = np.array([cen for cen in SGData['SGCen']],order='F')
2160   
2161    import pyspg
2162    Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
2163    h,k,l,f = Uniq
2164    Uniq=np.array(list(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq])))
2165    phi = f[:Nuniq]
2166    return iabsnt,mulp,Uniq,phi
2167   
2168def checkSSLaue(HKL,SGData,SSGData):
2169    #Laue check here - Toss HKL if outside unique Laue part
2170    h,k,l,m = HKL
2171    if SGData['SGLaue'] == '2/m':
2172        if SGData['SGUniq'] == 'a':
2173            if 'a' in SSGData['modSymb'] and h == 0 and m < 0:
2174                return False
2175            elif 'b' in SSGData['modSymb'] and k == 0 and l ==0 and m < 0:
2176                return False
2177            else:
2178                return True
2179        elif SGData['SGUniq'] == 'b':
2180            if 'b' in SSGData['modSymb'] and k == 0 and m < 0:
2181                return False
2182            elif 'a' in SSGData['modSymb'] and h == 0 and l ==0 and m < 0:
2183                return False
2184            else:
2185                return True
2186        elif SGData['SGUniq'] == 'c':
2187            if 'g' in SSGData['modSymb'] and l == 0 and m < 0:
2188                return False
2189            elif 'a' in SSGData['modSymb'] and h == 0 and k ==0 and m < 0:
2190                return False
2191            else:
2192                return True
2193    elif SGData['SGLaue'] == 'mmm':
2194        if 'a' in SSGData['modSymb']:
2195            if h == 0 and m < 0:
2196                return False
2197            else:
2198                return True
2199        elif 'b' in SSGData['modSymb']:
2200            if k == 0 and m < 0:
2201                return False
2202            else:
2203                return True
2204        elif 'g' in SSGData['modSymb']:
2205            if l == 0 and m < 0:
2206                return False
2207            else:
2208                return True
2209    else:   #tetragonal, trigonal, hexagonal (& triclinic?)
2210        if l == 0 and m < 0:
2211            return False
2212        else:
2213            return True
2214       
2215def checkHKLextc(HKL,SGData):
2216    '''
2217    Checks if reflection extinct - does not check centering
2218
2219    :param HKL:  [h,k,l]
2220    :param SGData: space group data obtained from SpcGroup
2221    :returns: True if extinct; False if allowed
2222
2223    '''
2224    Ops = SGData['SGOps']
2225    OpM = np.array([op[0] for op in Ops])
2226    OpT = np.array([op[1] for op in Ops])
2227    HKLS = np.array([HKL,-HKL])     #Freidel's Law
2228    DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,3))
2229    PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2230    for dhkl,phkl in zip(DHKL,PHKL)[1:]:    #skip identity
2231        if dhkl.any():
2232            continue
2233        else:
2234            if phkl%1.:
2235                return True
2236    return False
2237
2238def checkMagextc(HKL,SGData):
2239    '''
2240    Checks if reflection magnetically extinct; does fullcheck (centering, too)
2241    uses algorthm from Gallego, et al., J. Appl. Cryst. 45, 1236-1247 (2012)
2242
2243    :param HKL:  [h,k,l]
2244    :param SGData: space group data obtained from SpcGroup; must have magnetic symmetry SpnFlp data
2245    :returns: True if magnetically extinct; False if allowed (to match GenHKLf)
2246
2247    '''
2248    Ops = SGData['SGOps']
2249    Ncen = len(SGData['SGCen'])
2250    OpM = np.array([op[0] for op in Ops])
2251    OpT = np.array([op[1] for op in Ops])
2252    if SGData['SGInv'] and not SGData['SGFixed']:
2253        OpM = np.vstack((OpM,-OpM))
2254        OpT = np.vstack((OpT,-OpT))%1.
2255    OpM = np.reshape(np.array(list(OpM)*Ncen),(-1,3,3))
2256    OpT = np.reshape(np.array([OpT+cen for cen in SGData['SGCen']]),(-1,3))
2257    Spn = SGData['SpnFlp'][:len(OpM)]
2258    Mag = np.array([nl.det(opm) for opm in OpM])*Spn
2259    DHKL = np.reshape(np.inner(HKL,OpM),(-1,3))
2260    PHKL = np.reshape(np.cos(2.0*np.pi*np.inner(HKL,OpT))*Mag,(-1,))[:,nxs,nxs]*OpM     #compute T(R,theta) eq(7)
2261    Ftest = np.random.rand(3)       #random magnetic moment
2262    Psum = np.zeros(3)
2263    nsum = 0.
2264    nA = 0
2265    for dhkl,phkl in zip(DHKL,PHKL):
2266        if not np.allclose(dhkl,HKL):           #test for eq(5)
2267            continue
2268        else:
2269            nA += 1
2270            nsum += np.trace(phkl)          #eq(8)
2271            pterm = np.inner(Ftest,phkl)    #eq(9)
2272            Psum += pterm
2273    if nsum/nA > 1.:        #only need to look at nA=1 frok eq(8)
2274        return False
2275    if np.allclose(Psum,np.zeros(3)):
2276        return True
2277    else:
2278        if np.inner(HKL,Psum):
2279            return True
2280        return False
2281   
2282def checkSSextc(HKL,SSGData):
2283    Ops = SSGData['SSGOps']
2284    OpM = np.array([op[0] for op in Ops])
2285    OpT = np.array([op[1] for op in Ops])
2286    HKLS = np.array([HKL,-HKL])     #Freidel's Law
2287    DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,4))
2288    PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2289    for dhkl,phkl in list(zip(DHKL,PHKL))[1:]:    #skip identity
2290        if dhkl.any():
2291            continue
2292        else:
2293            if phkl%1.:
2294                return False
2295    return True
2296   
2297################################################################################
2298#### Site symmetry tables
2299################################################################################
2300     
2301OprName = {
2302    '-6643':       ['-1',1],'6479' :    ['2(z)',2],'-6479':     ['m(z)',3],
2303    '6481' :     ['m(y)',4],'-6481':    ['2(y)',5],'6641' :     ['m(x)',6],
2304    '-6641':     ['2(x)',7],'6591' :  ['m(+-0)',8],'-6591':   ['2(+-0)',9],
2305    '6531' :  ['m(110)',10],'-6531': ['2(110)',11],'6537' :    ['4(z)',12],
2306    '-6537':   ['-4(z)',13],'975'  : ['3(111)',14],'6456' :       ['3',15],
2307    '-489' :  ['3(+--)',16],'483'  : ['3(-+-)',17],'-969' :  ['3(--+)',18],
2308    '819'  :  ['m(+0-)',19],'-819' : ['2(+0-)',20],'2431' :  ['m(0+-)',21],
2309    '-2431':  ['2(0+-)',22],'-657' :  ['m(xz)',23],'657'  :   ['2(xz)',24],
2310    '1943' :   ['-4(x)',25],'-1943':   ['4(x)',26],'-2429':   ['m(yz)',27],
2311    '2429' :   ['2(yz)',28],'639'  :  ['-4(y)',29],'-639' :    ['4(y)',30],
2312    '-6484':   ['2(010)',4],'6484' :  ['m(010)',5],'-6668':   ['2(100)',6],
2313    '6668' :   ['m(100)',7],'-6454': ['2(120)',18],'6454' :  ['m(120)',19],
2314    '-6638':  ['2(210)',20],'6638' : ['m(210)',21],   #search in SytSym ends at m(210)
2315    '2223' : ['3(+++)2',39],
2316    '6538' :   ['6(z)1',40],'-2169':['3(--+)2',41],'2151' : ['3(+--)2',42],
2317    '2205' :['-3(-+-)2',43],'-2205':[' (-+-)2',44],'489'  :['-3(+--)1',45],
2318    '801'  :   ['4(y)1',46],'1945' :  ['4(x)3',47],'-6585': ['-4(z)3 ',48],
2319    '6585' :   ['4(z)3',49],'6584' :  ['3(z)2',50],'6666' :  ['6(z)5 ',51],
2320    '6643' :       ['1',52],'-801' : ['-4(y)1',53],'-1945': ['-4(x)3 ',54],
2321    '-6666':  ['-6(z)5',55],'-6538': ['-6(z)1',56],'-2223':['-3(+++)2',57],
2322    '-975' :['-3(+++)1',58],'-6456': ['-3(z)1',59],'-483' :['-3(-+-)1',60],
2323    '969'  :['-3(--+)1',61],'-6584': ['-3(z)2',62],'2169' :['-3(--+)2',63],
2324    '-2151':['-3(+--)2',64],   }                               
2325
2326KNsym = {
2327    '0'         :'    1   ','1'         :'   -1   ','64'        :'    2(x)','32'        :'    m(x)',
2328    '97'        :'  2/m(x)','16'        :'    2(y)','8'         :'    m(y)','25'        :'  2/m(y)',
2329    '2'         :'    2(z)','4'         :'    m(z)','7'         :'  2/m(z)','134217728' :'   2(yz)',
2330    '67108864'  :'   m(yz)','201326593' :' 2/m(yz)','2097152'   :'  2(0+-)','1048576'   :'  m(0+-)',
2331    '3145729'   :'2/m(0+-)','8388608'   :'   2(xz)','4194304'   :'   m(xz)','12582913'  :' 2/m(xz)',
2332    '524288'    :'  2(+0-)','262144'    :'  m(+0-)','796433'    :'2/m(+0-)','1024'      :'   2(xy)',
2333    '512'       :'   m(xy)','1537'      :' 2/m(xy)','256'       :'  2(+-0)','128'       :'  m(+-0)',
2334    '385'       :'2/m(+-0)','76'        :'  mm2(x)','52'        :'  mm2(y)','42'        :'  mm2(z)',
2335    '135266336' :' mm2(yz)','69206048'  :'mm2(0+-)','8650760'   :' mm2(xz)','4718600'   :'mm2(+0-)',
2336    '1156'      :' mm2(xy)','772'       :'mm2(+-0)','82'        :'  222   ','136314944' :'  222(x)',
2337    '8912912'   :'  222(y)','1282'      :'  222(z)','127'       :'  mmm   ','204472417' :'  mmm(x)',
2338    '13369369'  :'  mmm(y)','1927'      :'  mmm(z)','33554496'  :'  4(x)','16777280'  :' -4(x)',
2339    '50331745'  :'4/m(x)'  ,'169869394' :'422(x)','84934738'  :'-42m(x)','101711948' :'4mm(x)',
2340    '254804095' :'4/mmm(x)','536870928 ':'  4(y)','268435472' :' -4(y)','805306393' :'4/m(y)',
2341    '545783890' :'422(y)','272891986' :'-42m(y)','541327412' :'4mm(y)','818675839' :'4/mmm(y)',
2342    '2050'      :'  4(z)','4098'      :' -4(z)','6151'      :'4/m(z)','3410'      :'422(z)',
2343    '4818'      :'-42m(z)','2730'      :'4mm(z)','8191'      :'4/mmm(z)','8192'      :'  3(111)',
2344    '8193'      :' -3(111)','2629888'   :' 32(111)','1319040'   :' 3m(111)','3940737'   :'-3m(111)',
2345    '32768'     :'  3(+--)','32769'     :' -3(+--)','10519552'  :' 32(+--)','5276160'   :' 3m(+--)',
2346    '15762945'  :'-3m(+--)','65536'     :'  3(-+-)','65537'     :' -3(-+-)','134808576' :' 32(-+-)',
2347    '67437056'  :' 3m(-+-)','202180097' :'-3m(-+-)','131072'    :'  3(--+)','131073'    :' -3(--+)',
2348    '142737664' :' 32(--+)','71434368'  :' 3m(--+)','214040961' :'-3m(--+)','237650'    :'   23   ',
2349    '237695'    :'   m3   ','715894098' :'   432  ','358068946' :'  -43m  ','1073725439':'   m3m  ',
2350    '68157504'  :' mm2(d100)','4456464'   :' mm2(d010)','642'       :' mm2(d001)','153092172' :'-4m2(x)',
2351    '277348404' :'-4m2(y)','5418'      :'-4m2(z)','1075726335':'  6/mmm ','1074414420':'-6m2(100)',
2352    '1075070124':'-6m2(120)','1075069650':'   6mm  ','1074414890':'   622  ','1073758215':'   6/m  ',
2353    '1073758212':'   -6   ','1073758210':'    6   ','1073759865':'-3m(100)','1075724673':'-3m(120)',
2354    '1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
2355    '1073758209':'   -3   ','1073758208':'    3   ','1074135143':'mmm(100)','1075314719':'mmm(010)',
2356    '1073743751':'mmm(110)','1074004034':' mm2(z100)','1074790418':' mm2(z010)','1073742466':' mm2(z110)',
2357    '1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
2358    '1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
2359    '1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
2360    '1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
2361    '1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
2362    '1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
2363    '1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
2364    '1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':'   -1   ',
2365    }
2366
2367NXUPQsym = {
2368    '1'        :(28,29,28,28),'-1'       :( 1,29,28, 0),'2(x)'     :(12,18,12,25),'m(x)'     :(25,18,12,25),
2369    '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),
2370    '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),
2371    '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),
2372    '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),
2373    '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),
2374    '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),
2375    'mm2(x)'  :(12,10, 0,-1),'mm2(y)'   :(13,10, 0,-1),'mm2(z)'   :(14,10, 0,-1),
2376    'mm2(yz)'  :(10,13, 0,-1),'mm2(0+-)' :(11,13, 0,-1),'mm2(xz)'  :( 8,12, 0,-1),'mm2(+0-)' :( 9,12, 0,-1),
2377    'mm2(xy)'  :( 6,11, 0,-1),'222'      :( 1,10, 0,-1),'222(x)'   :( 1,13, 0,-1),
2378    '222(y)'   :( 1,12, 0,-1),'222(z)'   :( 1,11, 0,-1),'mmm'      :( 1,10, 0,-1),'mmm(x)'   :( 1,13, 0,-1),
2379    'mmm(y)'   :( 1,12, 0,-1),'mmm(z)'   :( 1,11, 0,-1),'4(x)'     :(12, 4,12, 0),'-4(x)'    :( 1, 4,12, 0),
2380    '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),
2381    '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),
2382    '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,),
2383    '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),
2384    '-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),
2385    '-3(111)'  :( 1, 5, 2, 0),'32(111)'  :( 1, 5, 0, 2),'3m(111)'  :( 2, 5, 0, 2),'-3m(111)' :( 1, 5, 0,-1),
2386    '3(+--)'   :( 5, 8, 5, 0),'-3(+--)'  :( 1, 8, 5, 0),'32(+--)'  :( 1, 8, 0, 5),'3m(+--)'  :( 5, 8, 0, 5),
2387    '-3m(+--)' :( 1, 8, 0,-1),'3(-+-)'   :( 4, 7, 4, 0),'-3(-+-)'  :( 1, 7, 4, 0),'32(-+-)'  :( 1, 7, 0, 4),
2388    '3m(-+-)'  :( 4, 7, 0, 4),'-3m(-+-)' :( 1, 7, 0,-1),'3(--+)'   :( 3, 6, 3, 0),'-3(--+)'  :( 1, 6, 3, 0),
2389    '32(--+)'  :( 1, 6, 0, 3),'3m(--+)'  :( 3, 6, 0, 3),'-3m(--+)' :( 1, 6, 0,-1),'23'       :( 1, 1, 0, 0),
2390    'm3'       :( 1, 1, 0, 0),'432'      :( 1, 1, 0, 0),'-43m'     :( 1, 1, 0, 0),'m3m'      :( 1, 1, 0, 0),
2391    'mm2(d100)':(12,13, 0,-1),'mm2(d010)':(13,12, 0,-1),'mm2(d001)':(14,11, 0,-1),'-4m2(x)'  :( 1, 4, 0,-1),
2392    '-4m2(y)'  :( 1, 3, 0,-1),'-4m2(z)'  :( 1, 2, 0,-1),'6/mmm'    :( 1, 9, 0,-1),'-6m2(100)':( 1, 9, 0,-1),
2393    '-6m2(120)':( 1, 9, 0,-1),'6mm'      :(14, 9, 0,-1),'622'      :( 1, 9, 0,-1),'6/m'      :( 1, 9,14,-1),
2394    '-6'       :( 1, 9,14, 0),'6'        :(14, 9,14, 0),'-3m(100)' :( 1, 9, 0,-1),'-3m(120)' :( 1, 9, 0,-1),
2395    '3m(100)'  :(14, 9, 0,14),'3m(120)'  :(14, 9, 0,14),'32(100)'  :( 1, 9, 0,14),'32(120)'  :( 1, 9, 0,14),
2396    '-3'       :( 1, 9,14, 0),'3'        :(14, 9,14, 0),'mmm(100)' :( 1,14, 0,-1),'mmm(010)' :( 1,15, 0,-1),
2397    'mmm(110)' :( 1,11, 0,-1),'mm2(z100)':(14,14, 0,-1),'mm2(z010)':(14,15, 0,-1),'mm2(z110)':(14,11, 0,-1),
2398    'mm2(100)' :(12,14, 0,-1),'mm2(010)' :(13,15, 0,-1),'mm2(110)' :( 6,11, 0,-1),'mm2(120)' :(15,14, 0,-1),
2399    'mm2(210)' :(16,15, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222(100)' :( 1,14, 0,-1),'222(010)' :( 1,15, 0,-1),
2400    '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),
2401    '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),
2402    'm(001)'   :(23,16,14,23),'m(100)'   :(26,25,12,26),'m(010)'   :(27,28,13,27),'m(110)'   :(18,19, 6,18),
2403    'm(120)'   :(24,27,15,24),'m(210)'   :(25,26,16,25),'m(+-0)'   :(17,20, 7,17),'2(001)'   :(14,16,14,23),
2404    '2(100)'   :(12,25,12,26),'2(010)'   :(13,28,13,27),'2(110)'   :( 6,19, 6,18),'2(120)'   :(15,27,15,24),
2405    '2(210)'   :(16,26,16,25),'2(+-0)'   :( 7,20, 7,17),'-1'       :( 1,29,28, 0)
2406    }
2407       
2408CSxinel = [[],      # 0th empty - indices are Fortran style
2409    [[0,0,0],[ 0.0, 0.0, 0.0]],      #1  0  0  0
2410    [[1,1,1],[ 1.0, 1.0, 1.0]],      #2  X  X  X
2411    [[1,1,1],[ 1.0, 1.0,-1.0]],      #3  X  X -X
2412    [[1,1,1],[ 1.0,-1.0, 1.0]],      #4  X -X  X
2413    [[1,1,1],[ 1.0,-1.0,-1.0]],      #5 -X  X  X
2414    [[1,1,0],[ 1.0, 1.0, 0.0]],      #6  X  X  0
2415    [[1,1,0],[ 1.0,-1.0, 0.0]],      #7  X -X  0
2416    [[1,0,1],[ 1.0, 0.0, 1.0]],      #8  X  0  X
2417    [[1,0,1],[ 1.0, 0.0,-1.0]],      #9  X  0 -X
2418    [[0,1,1],[ 0.0, 1.0, 1.0]],      #10  0  Y  Y
2419    [[0,1,1],[ 0.0, 1.0,-1.0]],      #11 0  Y -Y
2420    [[1,0,0],[ 1.0, 0.0, 0.0]],      #12  X  0  0
2421    [[0,1,0],[ 0.0, 1.0, 0.0]],      #13  0  Y  0
2422    [[0,0,1],[ 0.0, 0.0, 1.0]],      #14  0  0  Z
2423    [[1,1,0],[ 1.0, 2.0, 0.0]],      #15  X 2X  0
2424    [[1,1,0],[ 2.0, 1.0, 0.0]],      #16 2X  X  0
2425    [[1,1,2],[ 1.0, 1.0, 1.0]],      #17  X  X  Z
2426    [[1,1,2],[ 1.0,-1.0, 1.0]],      #18  X -X  Z
2427    [[1,2,1],[ 1.0, 1.0, 1.0]],      #19  X  Y  X
2428    [[1,2,1],[ 1.0, 1.0,-1.0]],      #20  X  Y -X
2429    [[1,2,2],[ 1.0, 1.0, 1.0]],      #21  X  Y  Y
2430    [[1,2,2],[ 1.0, 1.0,-1.0]],      #22  X  Y -Y
2431    [[1,2,0],[ 1.0, 1.0, 0.0]],      #23  X  Y  0
2432    [[1,0,2],[ 1.0, 0.0, 1.0]],      #24  X  0  Z
2433    [[0,1,2],[ 0.0, 1.0, 1.0]],      #25  0  Y  Z
2434    [[1,1,2],[ 1.0, 2.0, 1.0]],      #26  X 2X  Z
2435    [[1,1,2],[ 2.0, 1.0, 1.0]],      #27 2X  X  Z
2436    [[1,2,3],[ 1.0, 1.0, 1.0]],      #28  X  Y  Z
2437    ]
2438
2439CSuinel = [[],      # 0th empty - indices are Fortran style
2440    [[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
2441    [[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
2442    [[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
2443    [[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
2444    [[1,1,1,2,2,2],[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],[1,0,0,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]],    #5  A  A  A  D  D  D
2445    [[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
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]],    #7  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]],    #8  A  A  A  D  D -D
2448    [[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
2449    [[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
2450    [[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
2451    [[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
2452    [[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
2453    [[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
2454    [[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
2455    [[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
2456    [[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
2457    [[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
2458    [[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
2459    [[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
2460    [[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
2461    [[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
2462    [[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
2463    [[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
2464    [[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
2465    [[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
2466    [[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
2467    [[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
2468    [[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
2469    ]
2470   
2471################################################################################
2472#### Site symmetry routines
2473################################################################################
2474   
2475def GetOprPtrName(key):
2476    'Needs a doc string'
2477    try:
2478        oprName = OprName[key][0]
2479    except KeyError:
2480        return key
2481    return oprName.replace('(','').replace(')','')
2482
2483def GetOprPtrNumber(key):
2484    'Needs a doc string'
2485    try:
2486        return OprName[key][1]
2487    except KeyError:
2488        return key
2489
2490def GetOprName(key):
2491    'Needs a doc string'
2492    return OprName[key][0]
2493
2494def GetKNsym(key):
2495    'Needs a doc string'
2496    try:
2497        return KNsym[key].strip()
2498    except KeyError:
2499        return 'sp'
2500
2501def GetNXUPQsym(siteSym):
2502    '''       
2503    The codes XUPQ are for lookup of symmetry constraints for position(X), thermal parm(U) & magnetic moments (P & Q)
2504    '''
2505    return NXUPQsym[siteSym]
2506
2507def GetCSxinel(siteSym): 
2508    "returns Xyz terms, multipliers, GUI flags"
2509    indx = GetNXUPQsym(siteSym.strip())
2510    return CSxinel[indx[0]]
2511   
2512def GetCSuinel(siteSym):
2513    "returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
2514    indx = GetNXUPQsym(siteSym.strip())
2515    return CSuinel[indx[1]]
2516   
2517def GetCSpqinel(SpnFlp,dupDir): 
2518    "returns Mxyz terms, multipliers, GUI flags"
2519    CSI = [[1,2,3],[1.0,1.0,1.0]]
2520    for sopr in dupDir:
2521#        print (sopr,dupDir[sopr])
2522        opr = sopr.replace("'",'')
2523        indx = GetNXUPQsym(opr)
2524        if SpnFlp[dupDir[sopr]] > 0:
2525            csi = CSxinel[indx[2]]  #P
2526        else:
2527            csi = CSxinel[indx[3]]  #Q
2528#        print(opr,indx,csi,CSI)
2529        if not len(csi):
2530            return [[0,0,0],[0.,0.,0.]]
2531        for kcs in [0,1,2]:
2532            if csi[0][kcs] == 0 and CSI[0][kcs] != 0:
2533                jcs = CSI[0][kcs]
2534                for ics in [0,1,2]:
2535                    if CSI[0][ics] == jcs:
2536                        CSI[0][ics] = 0
2537                        CSI[1][ics] = 0.
2538                    elif CSI[0][ics] > jcs:
2539                        CSI[0][ics] = CSI[0][ics]-1
2540            elif (CSI[0][kcs] == csi[0][kcs]) and (CSI[1][kcs] != csi[1][kcs]):
2541                CSI[1][kcs] = csi[1][kcs]
2542            elif CSI[0][kcs] >= csi[0][kcs]:
2543                CSI[0][kcs] = min(CSI[0][kcs],csi[0][kcs])
2544                if CSI[1][kcs] != csi[1][kcs]:
2545                    if CSI[1][kcs] == 1.:
2546                        CSI[1][kcs] = csi[1][kcs]
2547#        print(CSI)
2548    return CSI
2549   
2550def getTauT(tau,sop,ssop,XYZ,wave=np.zeros(3)):
2551    phase = np.sum(XYZ*wave)
2552    ssopinv = nl.inv(ssop[0])
2553    mst = ssopinv[3][:3]
2554    epsinv = ssopinv[3][3]
2555    sdet = nl.det(sop[0])
2556    ssdet = nl.det(ssop[0])
2557    dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
2558    dT = 1.0
2559    if np.any(dtau%.5):
2560        sumdtau = np.sum(dtau%.5)
2561        dT = 0.
2562        if np.abs(sumdtau-.5) > 1.e-4:
2563            dT = np.tan(np.pi*sumdtau)
2564    tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]+phase)
2565    return sdet,ssdet,dtau,dT,tauT
2566   
2567def OpsfromStringOps(A,SGData,SSGData):
2568    SGOps = SGData['SGOps']
2569    SSGOps = SSGData['SSGOps']
2570    Ax = A.split('+')
2571    Ax[0] = int(Ax[0])
2572    iC = 1
2573    if Ax[0] < 0:
2574        iC = -1
2575    iAx = abs(Ax[0])
2576    nA = iAx%100-1
2577    nC = iAx//100
2578    unit = [0,0,0]
2579    if len(Ax) > 1:
2580        unit = eval('['+Ax[1]+']')
2581    return SGOps[nA],SSGOps[nA],iC,SGData['SGCen'][nC],unit
2582   
2583def GetSSfxuinel(waveType,Stype,nH,XYZ,SGData,SSGData,debug=False):
2584   
2585    def orderParms(CSI):
2586        parms = [0,]
2587        for csi in CSI:
2588            for i in [0,1,2]:
2589                if csi[i] not in parms:
2590                    parms.append(csi[i])
2591        for csi in CSI:
2592            for i in [0,1,2]:
2593                csi[i] = parms.index(csi[i])
2594        return CSI
2595       
2596    def fracCrenel(tau,Toff,Twid):
2597        Tau = (tau-Toff[:,nxs])%1.
2598        A = np.where(Tau<Twid[:,nxs],1.,0.)
2599        return A
2600       
2601    def fracFourier(tau,nH,fsin,fcos):
2602        SA = np.sin(2.*nH*np.pi*tau)
2603        CB = np.cos(2.*nH*np.pi*tau)
2604        A = SA[nxs,nxs,:]*fsin[:,:,nxs]
2605        B = CB[nxs,nxs,:]*fcos[:,:,nxs]
2606        return A+B
2607       
2608    def posFourier(tau,nH,psin,pcos):
2609        SA = np.sin(2*nH*np.pi*tau)
2610        CB = np.cos(2*nH*np.pi*tau)
2611        A = SA[nxs,nxs,:]*psin[:,:,nxs]
2612        B = CB[nxs,nxs,:]*pcos[:,:,nxs]
2613        return A+B   
2614
2615    def posZigZag(tau,Tmm,XYZmax):
2616        DT = Tmm[1]-Tmm[0]
2617        slopeUp = 2.*XYZmax/DT
2618        slopeDn = 2.*XYZmax/(1.-DT)
2619        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])
2620        return A
2621
2622    def posBlock(tau,Tmm,XYZmax):
2623        A = np.array([np.where(Tmm[0] < t <= Tmm[1],XYZmax,-XYZmax) for t in tau])
2624        return A
2625       
2626    def DoFrac():
2627        delt2 = np.eye(2)*0.001
2628        dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2629        dFTP = []
2630        if siteSym == '1':
2631            CSI = [[1,0],[2,0]],2*[[1.,0.],]
2632        elif siteSym == '-1':
2633            CSI = [[1,0],[0,0]],2*[[1.,0.],]
2634        else:
2635            FSC = np.ones(2,dtype='i')
2636            CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
2637            if 'Crenel' in waveType:
2638                dF = np.zeros_like(tau)
2639            else:
2640                dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2641            dFT = np.zeros_like(dF)
2642            dFTP = []
2643            for i in SdIndx:
2644                sop = Sop[i]
2645                ssop = SSop[i]           
2646                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2647                fsc = np.ones(2,dtype='i')
2648                if 'Crenel' in waveType:
2649                    dFT = np.zeros_like(tau)
2650                    fsc = [1,1]
2651                else:   #Fourier
2652                    dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
2653                    dFT = nl.det(sop[0])*dFT
2654                    dFT = dFT[:,np.argsort(tauT)]
2655                    dFT[0] *= ssdet
2656                    dFT[1] *= sdet
2657                    dFTP.append(dFT)
2658               
2659                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2660                        fsc = [1,1]
2661                        if dT:
2662                            CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
2663                        else:
2664                            CSI = [[[1,0],[0,0]],[[1.,0.],[0.,0.]]]
2665                        FSC = np.zeros(2,dtype='i')
2666                        return CSI,dF,dFTP
2667                    else:
2668                        for i in range(2):
2669                            if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
2670                                fsc[i] = 1
2671                            else:
2672                                fsc[i] = 0
2673                        FSC &= fsc
2674                        if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc)
2675            n = -1
2676            for i,F in enumerate(FSC):
2677                if F:
2678                    n += 1
2679                    CSI[0][i] = n+1
2680                    CSI[1][i] = 1.0
2681           
2682        return CSI,dF,dFTP
2683       
2684    def DoXYZ():
2685        delt5 = np.ones(5)*0.001
2686        delt6 = np.eye(6)*0.001
2687        if 'Fourier' in waveType:
2688            dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2689              #3x6x12 modulated position array (X,Spos,tau)& force positive
2690        elif waveType in ['ZigZag','Block']:
2691            if waveType == 'ZigZag':
2692                dX = posZigZag(tau,delt5[:2],delt5[2:])
2693            else:
2694                dX = posBlock(tau,delt5[:2],delt5[2:])
2695        dXTP = []
2696        if siteSym == '1':
2697            CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2698        elif siteSym == '-1':
2699            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.],]
2700        else:
2701            if 'Fourier' in waveType:
2702                CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2703            elif waveType in ['ZigZag','Block']:
2704                CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2705                    np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2706            XSC = np.ones(6,dtype='i')
2707            dXTP = []
2708            for i in SdIndx:
2709                sop = Sop[i]
2710                ssop = SSop[i]
2711                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2712                xsc = np.ones(6,dtype='i')
2713                if 'Fourier' in waveType:
2714                    dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,nxs,nxs]
2715                elif waveType == 'ZigZag':
2716                    dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2717                elif waveType == 'Block':
2718                    dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2719                dXT = np.inner(sop[0],dXT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2720                dXT = np.swapaxes(dXT,1,2)      # back to array(3x6x49)
2721                dXT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2722                dXTP.append(dXT)
2723                if waveType == 'Fourier':
2724                    for i in range(3):
2725                        if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2726                            xsc[i] = 0
2727                        if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2728                            xsc[i+3] = 0
2729                    if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2730                        xsc[3:6] = 0
2731                        CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2732                            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2733                        if dT:
2734                            if '(x)' in siteSym:
2735                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2736                                if 'm' in siteSym and len(SdIndx) == 1:
2737                                    CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2738                            elif '(y)' in siteSym:
2739                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2740                                if 'm' in siteSym and len(SdIndx) == 1:
2741                                    CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2742                            elif '(z)' in siteSym:
2743                                CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2744                                if 'm' in siteSym and len(SdIndx) == 1:
2745                                    CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2746                        else:
2747                            CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2748                    if '4/mmm' in laue:
2749                        if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2750                            if '(xy)' in siteSym:
2751                                CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2752                                if dT:
2753                                    CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2754                                else:
2755                                    CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2756                        if '(xy)' in siteSym or '(+-0)' in siteSym:
2757                            mul = 1
2758                            if '(+-0)' in siteSym:
2759                                mul = -1
2760                            if np.allclose(dX[0,0,:],dXT[1,0,:]):
2761                                CSI[0][3:5] = [[11,0,0],[11,0,0]]
2762                                CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2763                                xsc[3:5] = 0
2764                            if np.allclose(dX[0,3,:],dXT[0,4,:]):
2765                                CSI[0][:2] = [[12,0,0],[12,0,0]]
2766                                CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2767                                xsc[:2] = 0
2768                else:
2769                    for i in range(3):
2770                        if not np.allclose(dX[:,i],dXT[i,:,i]):
2771                            xsc[i] = 0
2772                XSC &= xsc
2773                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2774            if waveType == 'Fourier':
2775                n = -1
2776                if debug: print (XSC)
2777                for i,X in enumerate(XSC):
2778                    if X:
2779                        n += 1
2780                        CSI[0][i][0] = n+1
2781                        CSI[1][i][0] = 1.0
2782           
2783        return list(CSI),dX,dXTP
2784       
2785    def DoUij():
2786        delt12 = np.eye(12)*0.0001
2787        dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
2788        dUTP = []
2789        if siteSym == '1':
2790            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2791                [7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2792        elif siteSym == '-1':
2793            CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],   \
2794                6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2795        else:
2796            CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2797            USC = np.ones(12,dtype='i')
2798            dUTP = []
2799            dtau = 0.
2800            for i in SdIndx:
2801                sop = Sop[i]
2802                ssop = SSop[i]
2803                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2804                usc = np.ones(12,dtype='i')
2805                dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
2806                dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3)    #convert dUT to 12x49x3x3
2807                dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2808                dUT = np.array(U2Uij(dUijT))    #convert to 6x12x49
2809                dUT = dUT[:,:,np.argsort(tauT)]
2810                dUT[:,:6,:] *=(ssdet*sdet)
2811                dUTP.append(dUT)
2812                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2813                    if dT:
2814                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2815                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2816                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2817                        [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2818                    else:
2819                        CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 
2820                        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2821                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2822                        [0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2823                    if 'mm2(x)' in siteSym and dT:
2824                        CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2825                        USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2826                    elif '(xy)' in siteSym and dT:
2827                        CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2828                            [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2829                        CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2830                        USC = [1,1,1,1,1,1,1,1,1,1,1,1]                             
2831                    elif '(x)' in siteSym and dT:
2832                        CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2833                    elif '(y)' in siteSym and dT:
2834                        CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2835                    elif '(z)' in siteSym and dT:
2836                        CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2837                    for i in range(6):
2838                        if not USC[i]:
2839                            CSI[0][i] = [0,0,0]
2840                            CSI[1][i] = [0.,0.,0.]
2841                            CSI[0][i+6] = [0,0,0]
2842                            CSI[1][i+6] = [0.,0.,0.]
2843                else:                       
2844                    for i in range(6):
2845                        if not np.allclose(dU[i,i,:],dUT[i,i,:]):  #sin part
2846                            usc[i] = 0
2847                        if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]):   #cos part
2848                            usc[i+6] = 0
2849                    if np.any(dUT[1,0,:]):
2850                        if '4/m' in siteSym:
2851                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2852                            if ssop[1][3]:
2853                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2854                                usc[9] = 1
2855                            else:
2856                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2857                                usc[9] = 0
2858                        elif '4' in siteSym:
2859                            CSI[0][6:8] = [[12,0,0],[12,0,0]]
2860                            CSI[0][:2] = [[11,0,0],[11,0,0]]
2861                            if ssop[1][3]:
2862                                CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2863                                CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2864                                usc[2] = 0
2865                                usc[8] = 0
2866                                usc[3] = 1
2867                                usc[9] = 1
2868                            else:
2869                                CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2870                                CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2871                                usc[2] = 1
2872                                usc[8] = 1
2873                                usc[3] = 0               
2874                                usc[9] = 0
2875                        elif 'xy' in siteSym or '+-0' in siteSym:
2876                            if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2877                                CSI[0][4:6] = [[12,0,0],[12,0,0]]
2878                                CSI[0][6:8] = [[11,0,0],[11,0,0]]
2879                                CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2880                                CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2881                                usc[4:6] = 0
2882                                usc[6:8] = 0
2883                           
2884                    if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2885                USC &= usc
2886            if debug: print (USC)
2887            if not np.any(dtau%.5):
2888                n = -1
2889                for i,U in enumerate(USC):
2890                    if U:
2891                        n += 1
2892                        CSI[0][i][0] = n+1
2893                        CSI[1][i][0] = 1.0
2894   
2895        return list(CSI),dU,dUTP
2896   
2897    def DoMag():
2898        delt6 = np.eye(6)*0.001
2899        dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,nxs,nxs]
2900        dMTP = []
2901        CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2902        if siteSym == '1':
2903            CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2904        elif siteSym in ['-1','mmm',]:
2905            CSI = 3*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3*[[0.,0.,0.],]+3*[[1.,0.,0.],]
2906        elif siteSym in ['4(z)','422(z)']:
2907            CSI[0][0][0] = CSI[0][4][1] = 1
2908            CSI[1][0][0] = 1.0
2909            CSI[1][4][1] = -1.0
2910        elif siteSym in ['-4m2(z)','422(z)',]:
2911            CSI[0][5][0] = 1
2912            CSI[1][5][0] = 1.0
2913        elif siteSym in ['-32(100)','-3',]:
2914            CSI[0][2][0] = 1
2915            CSI[1][2][0] = 1.0
2916        elif siteSym in ['3',]:
2917            CSI[0][0][0] = CSI[0][3][0] = CSI[0][4][0] = 1
2918            CSI[1][0][0] = -np.sqrt(3.0)
2919            CSI[1][3][0] = 2.0
2920            CSI[1][4][0] = 1.0
2921        elif siteSym in ['622','2(100)','32(100)',]:
2922            CSI[0][0][0] = CSI[0][1][0] = CSI[0][3][0] = 1
2923            CSI[1][0][0] = 1.0
2924            CSI[1][1][0] = 2.0
2925            CSI[1][3][0] = np.sqrt(3.0)
2926        else:
2927              #3x6x12 modulated moment array (M,Spos,tau)& force positive
2928            CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2929            MSC = np.ones(6,dtype='i')
2930            dMTP = []
2931            for i in SdIndx:
2932                sop = Sop[i]
2933                ssop = SSop[i]
2934                sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2935                msc = np.ones(6,dtype='i')
2936                dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,nxs,nxs]
2937                dMT = np.inner(sop[0],dMT.T)    # X modulations array(3x6x49) -> array(3x49x6)
2938                dMT = np.swapaxes(dMT,1,2)      # back to array(3x6x49)
2939                dMT[:,:3,:] *= (ssdet*sdet)            # modify the sin component
2940                dMTP.append(dMT)
2941                for i in range(3):
2942                    if not np.allclose(dM[i,i,:],sdet*dMT[i,i,:]):
2943                        msc[i] = 0
2944                    if not np.allclose(dM[i,i+3,:],sdet*dMT[i,i+3,:]):
2945                        msc[i+3] = 0
2946                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2947                    msc[3:6] = 0
2948                    CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2949                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2950                    if dT:
2951                        if '(x)' in siteSym:
2952                            CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2953                            if 'm' in siteSym and len(SdIndx) == 1:
2954                                CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2955                        elif '(y)' in siteSym:
2956                            CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2957                            if 'm' in siteSym and len(SdIndx) == 1:
2958                                CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2959                        elif '(z)' in siteSym:
2960                            CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2961                            if 'm' in siteSym and len(SdIndx) == 1:
2962                                CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2963                    else:
2964                        CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2965                if '4/mmm' in laue:
2966                    if siteSym in ['4/mmm(z)',]:
2967                        CSI = 3*[[0,0,0],]+[[0,0,0],[0,0,0],[1,0,0]],3*[[0.,0.,0.],]+3*[[1.,0.,0.],]
2968                    if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2969                        if '(xy)' in siteSym:
2970                            CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2971                            if dT:
2972                                CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2973                            else:
2974                                CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2975                    if '(xy)' in siteSym or '(+-0)' in siteSym:
2976                        mul = 1
2977                        if '(+-0)' in siteSym:
2978                            mul = -1
2979                        if np.allclose(dM[0,0,:],dMT[1,0,:]):
2980                            CSI[0][3:5] = [[11,0,0],[11,0,0]]
2981                            CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2982                            msc[3:5] = 0
2983                        if np.allclose(dM[0,3,:],dMT[0,4,:]):
2984                            CSI[0][:2] = [[12,0,0],[12,0,0]]
2985                            CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2986                            msc[:2] = 0
2987                MSC &= msc
2988                if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2989            n = -1
2990            if debug: print (MSC)
2991            for i,M in enumerate(MSC):
2992                if M:
2993                    n += 1
2994                    CSI[0][i][0] = n+1
2995                    CSI[1][i][0] = 1.0
2996
2997        return list(CSI),dM,dMTP
2998       
2999    if debug: print ('super space group: '+SSGData['SSpGrp'])
3000    xyz = np.array(XYZ)%1.
3001    SGOps = copy.deepcopy(SGData['SGOps'])
3002    laue = SGData['SGLaue']
3003    siteSym = SytSym(XYZ,SGData)[0].strip()
3004    if debug: print ('siteSym: '+siteSym)
3005    SSGOps = copy.deepcopy(SSGData['SSGOps'])
3006    #expand ops to include inversions if any
3007    if SGData['SGInv'] and not SGData['SGFixed']:
3008        for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
3009            SGOps.append([-op[0],-op[1]%1.])
3010            SSGOps.append([-sop[0],-sop[1]%1.])
3011    #build set of sym ops around special position       
3012    SSop = []
3013    Sop = []
3014    Sdtau = []
3015    for iop,Op in enumerate(SGOps):         
3016        nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
3017        if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
3018            SSop.append(SSGOps[iop])
3019            Sop.append(SGOps[iop])
3020            ssopinv = nl.inv(SSGOps[iop][0])
3021            mst = ssopinv[3][:3]
3022            epsinv = ssopinv[3][3]
3023            Sdtau.append(np.sum(mst*(XYZ-SGOps[iop][1])-epsinv*SSGOps[iop][1][3]))
3024    SdIndx = np.argsort(np.array(Sdtau))     # just to do in sensible order
3025    if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
3026    #setup displacement arrays
3027    tau = np.linspace(-1,1,49,True)
3028    #make modulation arrays - one parameter at a time
3029    if Stype == 'Sfrac':
3030        CSI,dF,dFTP = DoFrac()
3031    elif Stype == 'Spos':
3032        CSI,dF,dFTP = DoXYZ()
3033        CSI[0] = orderParms(CSI[0])
3034    elif Stype == 'Sadp':
3035        CSI,dF,dFTP = DoUij()
3036        CSI[0] = orderParms(CSI[0]) 
3037    elif Stype == 'Smag':
3038        CSI,dF,dFTP = DoMag()
3039
3040    if debug:
3041        return CSI,dF,dFTP
3042    else:
3043        return CSI,[],[]
3044   
3045def MustrainNames(SGData):
3046    'Needs a doc string'
3047    laue = SGData['SGLaue']
3048    uniq = SGData['SGUniq']
3049    if laue in ['m3','m3m']:
3050        return ['S400','S220']
3051    elif laue in ['6/m','6/mmm','3m1','3']:
3052        return ['S400','S004','S202']
3053    elif laue in ['31m',]:
3054        return ['S400','S004','S202','S301']
3055    elif laue in ['3R','3mR']:
3056        return ['S400','S220','S310','S211']
3057    elif laue in ['4/m','4/mmm']:
3058        return ['S400','S004','S220','S022']
3059    elif laue in ['mmm']:
3060        return ['S400','S040','S004','S220','S202','S022']
3061    elif laue in ['2/m']:
3062        SHKL = ['S400','S040','S004','S220','S202','S022']
3063        if uniq == 'a':
3064            SHKL += ['S013','S031','S211']
3065        elif uniq == 'b':
3066            SHKL += ['S301','S103','S121']
3067        elif uniq == 'c':
3068            SHKL += ['S130','S310','S112']
3069        return SHKL
3070    else:
3071        SHKL = ['S400','S040','S004','S220','S202','S022']
3072        SHKL += ['S310','S103','S031','S130','S301','S013']
3073        SHKL += ['S211','S121','S112']
3074        return SHKL
3075       
3076def HStrainVals(HSvals,SGData):
3077    laue = SGData['SGLaue']
3078    uniq = SGData['SGUniq']
3079    DIJ = np.zeros(6)
3080    if laue in ['m3','m3m']:
3081        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
3082    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3083        DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
3084    elif laue in ['3R','3mR']:
3085        DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
3086    elif laue in ['4/m','4/mmm']:
3087        DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
3088    elif laue in ['mmm']:
3089        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3090    elif laue in ['2/m']:
3091        DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3092        if uniq == 'a':
3093            DIJ[5] = HSvals[3]
3094        elif uniq == 'b':
3095            DIJ[4] = HSvals[3]
3096        elif uniq == 'c':
3097            DIJ[3] = HSvals[3]
3098    else:
3099        DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
3100    return DIJ
3101
3102def HStrainNames(SGData):
3103    'Needs a doc string'
3104    laue = SGData['SGLaue']
3105    uniq = SGData['SGUniq']
3106    if laue in ['m3','m3m']:
3107        return ['D11','eA']         #add cubic strain term
3108    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3109        return ['D11','D33']
3110    elif laue in ['3R','3mR']:
3111        return ['D11','D12']
3112    elif laue in ['4/m','4/mmm']:
3113        return ['D11','D33']
3114    elif laue in ['mmm']:
3115        return ['D11','D22','D33']
3116    elif laue in ['2/m']:
3117        Dij = ['D11','D22','D33']
3118        if uniq == 'a':
3119            Dij += ['D23']
3120        elif uniq == 'b':
3121            Dij += ['D13']
3122        elif uniq == 'c':
3123            Dij += ['D12']
3124        return Dij
3125    else:
3126        Dij = ['D11','D22','D33','D12','D13','D23']
3127        return Dij
3128   
3129def MustrainCoeff(HKL,SGData):
3130    'Needs a doc string'
3131    #NB: order of terms is the same as returned by MustrainNames
3132    laue = SGData['SGLaue']
3133    uniq = SGData['SGUniq']
3134    h,k,l = HKL
3135    Strm = []
3136    if laue in ['m3','m3m']:
3137        Strm.append(h**4+k**4+l**4)
3138        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
3139    elif laue in ['6/m','6/mmm','3m1','3']:
3140        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
3141        Strm.append(l**4)
3142        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
3143    elif laue in ['31m',]:
3144        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
3145        Strm.append(l**4)
3146        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
3147        Strm.append(4.0*l*h**3)
3148    elif laue in ['3R','3mR']:
3149        Strm.append(h**4+k**4+l**4)
3150        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
3151        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
3152        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
3153    elif laue in ['4/m','4/mmm']:
3154        Strm.append(h**4+k**4)
3155        Strm.append(l**4)
3156        Strm.append(3.0*(h*k)**2)
3157        Strm.append(3.0*((h*l)**2+(k*l)**2))
3158    elif laue in ['mmm']:
3159        Strm.append(h**4)
3160        Strm.append(k**4)
3161        Strm.append(l**4)
3162        Strm.append(3.0*(h*k)**2)
3163        Strm.append(3.0*(h*l)**2)
3164        Strm.append(3.0*(k*l)**2)
3165    elif laue in ['2/m']:
3166        Strm.append(h**4)
3167        Strm.append(k**4)
3168        Strm.append(l**4)
3169        Strm.append(3.0*(h*k)**2)
3170        Strm.append(3.0*(h*l)**2)
3171        Strm.append(3.0*(k*l)**2)
3172        if uniq == 'a':
3173            Strm.append(2.0*k*l**3)
3174            Strm.append(2.0*l*k**3)
3175            Strm.append(4.0*k*l*h**2)
3176        elif uniq == 'b':
3177            Strm.append(2.0*l*h**3)
3178            Strm.append(2.0*h*l**3)
3179            Strm.append(4.0*h*l*k**2)
3180        elif uniq == 'c':
3181            Strm.append(2.0*h*k**3)
3182            Strm.append(2.0*k*h**3)
3183            Strm.append(4.0*h*k*l**2)
3184    else:
3185        Strm.append(h**4)
3186        Strm.append(k**4)
3187        Strm.append(l**4)
3188        Strm.append(3.0*(h*k)**2)
3189        Strm.append(3.0*(h*l)**2)
3190        Strm.append(3.0*(k*l)**2)
3191        Strm.append(2.0*k*h**3)
3192        Strm.append(2.0*h*l**3)
3193        Strm.append(2.0*l*k**3)
3194        Strm.append(2.0*h*k**3)
3195        Strm.append(2.0*l*h**3)
3196        Strm.append(2.0*k*l**3)
3197        Strm.append(4.0*k*l*h**2)
3198        Strm.append(4.0*h*l*k**2)
3199        Strm.append(4.0*k*h*l**2)
3200    return Strm
3201
3202def MuShklMean(SGData,Amat,Shkl):
3203   
3204    def genMustrain(xyz,Shkl):
3205        uvw = np.inner(Amat.T,xyz)
3206        Strm = np.array(MustrainCoeff(uvw,SGData))
3207        Sum = np.sum(np.multiply(Shkl,Strm))
3208        Sum = np.where(Sum > 0.01,Sum,0.01)
3209        Sum = np.sqrt(Sum)
3210        return Sum*xyz
3211       
3212    PHI = np.linspace(0.,360.,30,True)
3213    PSI = np.linspace(0.,180.,30,True)
3214    X = np.outer(npcosd(PHI),npsind(PSI))
3215    Y = np.outer(npsind(PHI),npsind(PSI))
3216    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3217    XYZ = np.dstack((X,Y,Z))
3218    XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
3219    return np.sqrt(np.sum(XYZ**2)/900.)
3220   
3221def Muiso2Shkl(muiso,SGData,cell):
3222    "this is to convert isotropic mustrain to generalized Shkls"
3223    import GSASIIlattice as G2lat
3224    A = G2lat.cell2AB(cell)[0]
3225   
3226    def minMus(Shkl,muiso,H,SGData,A):
3227        U = np.inner(A.T,H)
3228        S = np.array(MustrainCoeff(U,SGData))
3229        nS = S.shape[0]
3230        Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:nS,nxs]),axis=0))
3231        rad = np.sqrt(np.sum((Sum[:,nxs]*H)**2,axis=1))
3232        return (muiso-rad)**2
3233       
3234    laue = SGData['SGLaue']
3235    PHI = np.linspace(0.,360.,60,True)
3236    PSI = np.linspace(0.,180.,60,True)
3237    X = np.outer(npsind(PHI),npsind(PSI))
3238    Y = np.outer(npcosd(PHI),npsind(PSI))
3239    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3240    HKL = np.dstack((X,Y,Z))
3241    if laue in ['m3','m3m']:
3242        S0 = [1000.,1000.]
3243    elif laue in ['6/m','6/mmm']:
3244        S0 = [1000.,1000.,1000.]
3245    elif laue in ['31m','3','3m1']:
3246        S0 = [1000.,1000.,1000.]
3247    elif laue in ['3R','3mR']:
3248        S0 = [1000.,1000.,1000.,1000.]
3249    elif laue in ['4/m','4/mmm']:
3250        S0 = [1000.,1000.,1000.,1000.]
3251    elif laue in ['mmm']:
3252        S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
3253    elif laue in ['2/m']:
3254        S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
3255    else:
3256        S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000., 
3257            1000.,1000.,0.,0.,0.]
3258    S0 = np.array(S0)
3259    HKL = np.reshape(HKL,(-1,3))
3260    result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
3261    return result[0]
3262       
3263def PackRot(SGOps):
3264    IRT = []
3265    for ops in SGOps:
3266        M = ops[0]
3267        irt = 0
3268        for j in range(2,-1,-1):
3269            for k in range(2,-1,-1):
3270                irt *= 3
3271                irt += M[k][j]
3272        IRT.append(int(irt))
3273    return IRT
3274       
3275def SytSym(XYZ,SGData):
3276    '''
3277    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
3278
3279    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
3280    :param SGData: from SpcGroup
3281    :Returns: a four element tuple:
3282
3283     * The 1st element is a code for the site symmetry (see GetKNsym)
3284     * The 2nd element is the site multiplicity
3285     * Ndup number of overlapping operators
3286     * dupDir Dict - dictionary of overlapping operators
3287
3288    '''
3289    if SGData['SpGrp'] == 'P 1':
3290        return '1',1,1,{}
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(np.array(XYZ)%1.,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    try:
3330        Mult = len(SGData['SGOps'])*Ncen*inv//Jdup
3331    except: # patch because Jdup is not getting incremented for most atoms!
3332        Mult = 0
3333       
3334    return GetKNsym(str(Isym)),Mult,Ndup,dupDir
3335   
3336def MagSytSym(SytSym,dupDir,SGData):
3337    '''
3338    site sym operations: 1,-1,2,3,-3,4,-4,6,-6,m need to be marked if spin inversion
3339    '''
3340    SGData['GenSym'],SGData['GenFlg'] = GetGenSym(SGData)[:2]
3341#    print('SGPtGrp',SGData['SGPtGrp'],'SytSym',SytSym,'MagSpGrp',SGData['MagSpGrp'])
3342#    print('dupDir',dupDir)
3343    SplitSytSym = SytSym.split('(')
3344    if SGData['SGGray']:
3345        return SytSym+"1'"
3346    if SytSym in ['1','-1','2','3','6','m']:       #axial position
3347        return SytSym
3348    if SplitSytSym[0] == SGData['SGPtGrp']:     #simple cases
3349        try:
3350            MagSytSym = SGData['MagSpGrp'].split()[1]
3351        except IndexError:
3352            MagSytSym = SGData['MagSpGrp'][1:].strip("1'")
3353        if len(SplitSytSym) > 1:
3354            MagSytSym += '('+SplitSytSym[1]
3355        return MagSytSym
3356    if len(dupDir) == 1:
3357        return list(dupDir.keys())[0]
3358   
3359   
3360    if '2/m' in SytSym:         #done I think; last 2wo might be not needed
3361        ops = {'(x)':['2(x)','m(x)'],'(y)':['2(y)','m(y)'],'(z)':['2(z)','m(z)'],
3362               '(100)':['2(100)','m(100)'],'(010)':['2(010)','m(010)'],'(001)':['2(001)','m(001)'],
3363               '(120)':['2(120)','m(120)'],'(210)':['2(210)','m(210)'],'(+-0)':['2(+-0)','m(+-0)'],
3364               '(110)':['2(110)','m(110)']}
3365   
3366    elif '4/mmm' in SytSym:
3367        ops = {'(x)':['4(x)','m(x)','m(y)','m(0+-)'],   #m(0+-) for cubic m3m?
3368               '(y)':['4(y)','m(y)','m(z)','m(+0-)'],   #m(+0-)
3369               '(z)':['4(z)','m(z)','m(x)','m(+-0)']}   #m(+-0)
3370    elif '4mm' in SytSym:
3371        ops = {'(x)':['4(x)','m(y)','m(yz)'],'(y)':['4(y)','m(z)','m(xz)'],'(z)':['4(z)','m(x)','m(110)']}
3372    elif '422' in SytSym:
3373        ops = {'(x)':['4(x)','2(y)','2(yz)'],'(y)':['4(y)','2(z)','2(xz)'],'(z)':['4(z)','2(x)','2(110)']}
3374    elif '-4m2' in SytSym:
3375        ops = {'(x)':['-4(x)','m(x)','2(yz)'],'(y)':['-4(y)','m(y)','2(xz)'],'(z)':['-4(z)','m(z)','2(110)']}
3376    elif '-42m' in SytSym:
3377        ops = {'(x)':['-4(x)','2(y)','m(yz)'],'(y)':['-4(y)','2(z)','m(xz)'],'(z)':['-4(z)','2(x)','m(110)']}
3378    elif '-4' in SytSym:
3379        ops = {'(x)':['-4(x)',],'(y)':['-4(y)',],'(z)':['-4(z)',],}
3380    elif '4' in SytSym:
3381        ops = {'(x)':['4(x)',],'(y)':['4(y)',],'(z)':['4(z)',],}
3382
3383    elif '222' in SytSym:
3384        ops = {'':['2(x)','2(y)','2(z)'],
3385                   '(x)':['2(y)','2(z)','2(x)'],'(y)':['2(x)','2(z)','2(y)'],'(z)':['2(x)','2(y)','2(z)'],
3386                   '(100)':['2(z)','2(100)','2(120)',],'(010)':['2(z)','2(010)','2(210)',],
3387                   '(110)':['2(z)','2(110)','2(+-0)',],}
3388    elif 'mm2' in SytSym:
3389        ops = {'(x)':['m(y)','m(z)','2(x)'],'(y)':['m(x)','m(z)','2(y)'],'(z)':['m(x)','m(y)','2(z)'],
3390               '(xy)':['m(+-0)','m(z)','2(110)'],'(yz)':['m(0+-)','m(xz)','2(yz)'],     #not 2(xy)!
3391               '(xz)':['m(+0-)','m(y)','2(xz)'],'(z100)':['m(100)','m(120)','2(z)'],
3392               '(z010)':['m(010)','m(210)','2(z)'],'(z110)':['m(110)','m(+-0)','2(z)'],
3393               '(+-0)':[ 'm(110)','m(z)','2(+-0)'],'(d100)':['m(yz)','m(0+-)','2(xz)'],
3394               '(d010)':['m(xz)','m(+0-)','2(y)'],'(d001)':['m(110)','m(+-0)','2(z)'],
3395               '(210)':['m(z)','m(010)','2(210)'],'(120)':['m(z)','m(100)','2(120)'],
3396               '(100)':['m(z)','m(120)','2(100)',],'(010)':['m(z)','m(210)','2(010)',],
3397               '(110)':['m(z)','m(+-0)','2(110)',],}
3398    elif 'mmm' in SytSym:
3399        ops = {'':['m(x)','m(y)','m(z)'],
3400                   '(100)':['m(z)','m(100)','m(120)',],'(010)':['m(z)','m(010)','m(210)',],
3401                   '(110)':['m(z)','m(110)','m(+-0)',],
3402                   '(x)':['m(x)','m(y)','m(z)'],'(y)':['m(x)','m(y)','m(z)'],'(z)':['m(x)','m(y)','m(z)'],}
3403       
3404    elif '32' in SytSym:
3405        ops = {'(120)':['3','2(120)',],'(100)':['3','2(100)'],'(111)':['3(111)','2(x)']}
3406    elif '23' in SytSym:
3407        ops = {'':['2(x)','3(111)']}
3408    elif 'm3' in SytSym:
3409        ops = {'(100)':['(+-0)',],'(+--)':[],'(-+-)':[],'(--+)':[]}
3410    elif '3m' in SytSym:
3411        ops = {'(111)':['3(111)','m(+-0)',],'(+--)':['3(+--)','m(0+-)',],
3412               '(-+-)':['3(-+-)','m(+0-)',],'(--+)':['3(--+)','m(+-0)',],
3413               '(100)':['3','m(100)'],'(120)':['3','m(210)',]}
3414   
3415    if SytSym.split('(')[0] in ['6/m','6mm','-6m2','622','-6','-3','-3m','-43m',]:     #not simple cases
3416        MagSytSym = SytSym
3417        if "-1'" in dupDir:
3418            if '-6' in SytSym:
3419                MagSytSym = MagSytSym.replace('-6',"-6'")
3420            elif '-3m' in SytSym:
3421                MagSytSym = MagSytSym.replace('-3m',"-3'm'")
3422            elif '-3' in SytSym:
3423                MagSytSym = MagSytSym.replace('-3',"-3'")
3424        elif '-6m2' in SytSym:
3425            if "m'(110)" in dupDir:
3426                MagSytSym = "-6m'2'("+SytSym.split('(')[1]
3427        elif '6/m' in SytSym:
3428            if "m'(z)" in dupDir:
3429                MagSytSym = "6'/m'"
3430        elif '6mm' in SytSym:
3431            if "m'(110)" in dupDir:
3432                MagSytSym = "6'm'm"
3433        elif '-43m' in SytSym:
3434            if "m'(110)" in dupDir:
3435                MagSytSym = "-43m'"
3436        return MagSytSym
3437    try:
3438        axis = '('+SytSym.split('(')[1]
3439    except IndexError:
3440        axis = ''
3441    MagSytSym = ''
3442    for m in ops[axis]:
3443        if m in dupDir:
3444            MagSytSym += m.split('(')[0]
3445        else:
3446            MagSytSym += m.split('(')[0]+"'"
3447        if '2/m' in SytSym and '2' in m:
3448            MagSytSym += '/'
3449        if '-3/m' in SytSym:
3450            MagSytSym = '-'+MagSytSym
3451       
3452    MagSytSym += axis
3453# some exceptions & special rules         
3454    if MagSytSym == "4'/m'm'm'": MagSytSym = "4/m'm'm'"
3455    return MagSytSym
3456   
3457#    if len(GenSym) == 3:
3458#        if SGSpin[1] < 0:
3459#            if 'mm2' in SytSym:
3460#                MagSytSym = "m'm'2"+'('+SplitSytSym[1]
3461#            else:   #bad rule for I41/a
3462#                MagSytSym = SplitSytSym[0]+"'"
3463#                if len(SplitSytSym) > 1:
3464#                    MagSytSym += '('+SplitSytSym[1]
3465#        else:
3466#            MagSytSym = SytSym
3467#        if len(SplitSytSym) >1:
3468#            if "-4'"+'('+SplitSytSym[1] in dupDir:
3469#                MagSytSym = MagSytSym.replace('-4',"-4'")
3470#            if "-6'"+'('+SplitSytSym[1] in dupDir:
3471#                MagSytSym = MagSytSym.replace('-6',"-6'")
3472#        return MagSytSym
3473#           
3474    return SytSym
3475
3476def UpdateSytSym(Phase):
3477    ''' Update site symmetry/site multiplicity after space group/BNS lattice change
3478    '''
3479    generalData = Phase['General']
3480    SGData = generalData['SGData']
3481    Atoms = Phase['Atoms']
3482    cx,ct,cs,cia = generalData['AtomPtrs']
3483    for atom in Atoms:
3484        XYZ = atom[cx:cx+3]
3485        sytsym,Mult = SytSym(XYZ,SGData)[:2]
3486        sytSym,Mul,Nop,dupDir = SytSym(XYZ,SGData)
3487        atom[cs] = sytsym
3488        if generalData['Type'] == 'magnetic':
3489            magSytSym = MagSytSym(sytSym,dupDir,SGData)
3490            atom[cs] = magSytSym
3491        atom[cs+1] = Mult
3492    return
3493   
3494def ElemPosition(SGData):
3495    ''' Under development.
3496    Object here is to return a list of symmetry element types and locations suitable
3497    for say drawing them.
3498    So far I have the element type... getting all possible locations without lookup may be impossible!
3499    '''
3500    Inv = SGData['SGInv']
3501    eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
3502    # get operators & expand if centrosymmetric
3503    SymElements = []
3504    Ops = SGData['SGOps']
3505    opM = np.array([op[0].T for op in Ops])
3506    opT = np.array([op[1] for op in Ops])
3507    if Inv:
3508        opM = np.concatenate((opM,-opM))
3509        opT = np.concatenate((opT,-opT))
3510    opMT = list(zip(opM,opT))
3511    for M,T in opMT[1:]:        #skip I
3512        Dt = int(nl.det(M))
3513        Tr = int(np.trace(M))
3514        Dt = -(Dt-1)//2
3515        Es = eleSym[Tr][Dt]
3516        if Dt:              #rotation-inversion
3517            I = np.eye(3)
3518            if Tr == 1:     #mirrors/glides
3519                if np.any(T):       #glide
3520                    M2 = np.inner(M,M)
3521                    MT = np.inner(M,T)+T
3522                    print ('glide',Es,MT)
3523                    print (M2)
3524                else:               #mirror
3525                    print ('mirror',Es,T)
3526                    print (I-M)
3527                X = [-1,-1,-1]
3528            elif Tr == -3:  # pure inversion
3529                X = np.inner(nl.inv(I-M),T)
3530                print ('inversion',Es,X)
3531            else:           #other rotation-inversion
3532                M2 = np.inner(M,M)
3533                MT = np.inner(M,T)+T
3534                print ('rot-inv',Es,MT)
3535                print (M2)
3536                X = [-1,-1,-1]
3537        else:               #rotations
3538            print ('rotation',Es)
3539            X = [-1,-1,-1]
3540        SymElements.append([Es,X])
3541       
3542    return SymElements
3543   
3544def ApplyStringOps(A,SGData,X,Uij=[]):
3545    'Needs a doc string'
3546    SGOps = SGData['SGOps']
3547    SGCen = SGData['SGCen']
3548    Ax = A.split('+')
3549    Ax[0] = int(Ax[0])
3550    iC = 1
3551    if Ax[0] < 0:
3552        iC = -1
3553    Ax[0] = abs(Ax[0])
3554    nA = Ax[0]%100-1
3555    cA = Ax[0]//100
3556    Cen = SGCen[cA]
3557    M,T = SGOps[nA]
3558    if len(Ax)>1:
3559        cellA = Ax[1].split(',')
3560        cellA = np.array([int(a) for a in cellA])
3561    else:
3562        cellA = np.zeros(3)
3563    newX = Cen+iC*(np.inner(M,X).T+T)+cellA
3564    if len(Uij):
3565        U = Uij2U(Uij)
3566        U = np.inner(M,np.inner(U,M).T)
3567        newUij = U2Uij(U)
3568        return [newX,newUij]
3569    else:
3570        return newX
3571       
3572def ApplyStringOpsMom(A,SGData,SSGData,Mom):
3573    '''Applies string operations to modulated magnetic moment components used in drawing
3574    Drawing matches Bilbao MVISUALIZE
3575    '''
3576    SGOps = SGData['SGOps']
3577    Ax = A.split('+')
3578    Ax[0] = int(Ax[0])
3579    iAx = abs(Ax[0])
3580    nA = iAx%100-1
3581    if SGData['SGInv'] and not SGData['SGFixed']:
3582        nC = 2*len(SGOps)*(iAx//100)
3583    else:
3584        nC = len(SGOps)*(iAx//100)
3585    NA = nA
3586    if Ax[0] < 0:
3587        NA += len(SGOps)
3588    M,T = SGOps[nA]
3589    newMom = np.inner(Mom,M).T*SGData['SpnFlp'][NA+nC]
3590    if SSGData is not None:
3591        if SSGData['SSGCen'][iAx//100][3]:     #flip spin for BNS centered atoms
3592            newMom *= -1.
3593    return newMom
3594       
3595def StringOpsProd(A,B,SGData):
3596    """
3597    Find A*B where A & B are in strings '-' + '100*c+n' + '+ijk'
3598    where '-' indicates inversion, c(>0) is the cell centering operator,
3599    n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3600    Should return resultant string - C. SGData - dictionary using entries:
3601
3602       *  'SGCen': cell centering vectors [0,0,0] at least
3603       *  'SGOps': symmetry operations as [M,T] so that M*x+T = x'