source: trunk/GSASIIspc.py @ 3245

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

fix problem with incommensurate magnetic structures from CIFhklReader
make "Edit Body" as "Edit Residue Body" in RB menu
change protein validation bar plots to be vertical
change 'twin' to 'flag' in Reflection List table & add '-2' as 'free' (for proteins - not yet implemented)
fix PDB phase import to get atom type from 76:78 of ATOM record
change importer names for single crystal TOF data to be more explicit (SNS vs ISIS)
change cif structure factor importer - F2, sig(F2) & Fcalc now allowed

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