source: trunk/GSASIIspc.py @ 3293

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

Fix problems with GenAtom?

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