source: trunk/GSASIIspc.py @ 3340

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

make testDeriv py3 compatible
fix problem with Reload draw atoms - mag moments
add |Mag| to lst output of magnetic moments
correct errors in mag moment structure factor & derivatives

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