source: trunk/GSASIIspc.py @ 3234

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

new G2spc.SSChoice routine finds all allowed super space group symbols for a given space group.
Installed to drive ss symbol choices. Simplifies SSpcGroup.
Fix python 3 bug in covariance plot routine

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