source: trunk/GSASIIspc.py @ 3295

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

remove superfluous calls to SpcGroup? after SGData has been defined
add 3rd row magnetic form factors
add contour plots for projections of mustrain, size, % 3D pole distribution
improve quality of these surfaces

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