source: trunk/GSASIIspc.py @ 3297

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

progress on magnetism - BNS system ok for monoclinic, orthorhombic, tetragonal, hexagonal & trigonal
a few exceptions remain (cubic, rhombohedral, etc.)

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