source: trunk/GSASIIspc.py @ 3242

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

implement calc of mean mustrain from generalized mustrain coefficients

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