source: trunk/GSASIIspc.py @ 3300

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

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