source: trunk/GSASIIspc.py @ 3418

Last change on this file since 3418 was 3418, checked in by vondreele, 5 years ago

move SGMagSpinBox from G2phsGUI to G2ctrlGUI
add selection of spin choices in Unit Cells List for testing magnetic space groups with neutron data only - test TBD
fix bug in ApplyBNSlatt for non BNS centered lattice choice

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