source: trunk/GSASIIspc.py @ 3400

Last change on this file since 3400 was 3400, checked in by vondreele, 4 years ago

fix space group operator interpreter for magnetic mcif files (only place used!)

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