source: trunk/GSASIIspc.py @ 530

Last change on this file since 530 was 530, checked in by vondreele, 11 years ago

small fix to HessianLSQ about # of cycles
finish map peak search (no GUI/plot output yet)
  • Property svn:keywords set to Date Author Revision URL Id
File size: 52.0 KB
Line 
1"GSASII - Space group interpretion routines"
2
3########### SVN repository information ###################
4# $Date: 2012-04-06 18:30:45 +0000 (Fri, 06 Apr 2012) $
5# $Author: vondreele $
6# $Revision: 530 $
7# $URL: trunk/GSASIIspc.py $
8# $Id: GSASIIspc.py 530 2012-04-06 18:30:45Z vondreele $
9########### SVN repository information ###################
10import numpy as np
11import numpy.ma as ma
12import numpy.linalg as nl
13import math
14import sys
15import os.path as ospath
16
17import GSASIIpath
18import pyspg
19
20def SpcGroup(SGSymbol):
21    '''
22    Determines cell and symmetry information from a short H-M space group name
23    input:
24        SGSymbol - space group symbol (string) with spaces between axial fields
25    returns:
26        SGError = 0 for no errors; >0 for errors (see SGErrors below for details)
27        SGData - dictionary with entries:
28             'SpGrp': space group symbol slightly cleaned up
29             'Laue':  one of '-1','2/m','mmm','4/m','4/mmm','3R','3mR','3',
30                      '3m1','31m','6/m','6/mmm','m3','m3m'
31             'SGInv': boolean; True if centrosymmetric, False if not
32             'SGLatt': one of 'P','A','B','C','I','F','R'
33             'SGUniq': one of 'a','b','c' if monoclinic, '' otherwise
34             'SGCen': cell centering vectors [0,0,0] at least
35             'SGOps': symmetry operations as [M,T] so that M*x+T = x'
36             'SGSys': one of 'triclinic','monoclinic','orthorhombic','tetragonal','rhombohedral','trigonal','hexagonal','cubic'
37       '''
38    LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
39    LattSym = ('P','A','B','C','I','F','R')
40    UniqSym = ('','','a','b','c','',)
41    SysSym = ('triclinic','monoclinic','orthorhombic','tetragonal','rhombohedral','trigonal','hexagonal','cubic')
42    SGData = {}
43    SGData['SpGrp'] = SGSymbol.strip().lower().capitalize()
44    SGInfo = pyspg.sgforpy(SGSymbol)
45    SGData['SGLaue'] = LaueSym[SGInfo[0]-1]
46    SGData['SGInv'] = bool(SGInfo[1])
47    SGData['SGLatt'] = LattSym[SGInfo[2]-1]
48    SGData['SGUniq'] = UniqSym[SGInfo[3]+1]
49    if SGData['SGLatt'] == 'P':
50        SGData['SGCen'] = np.array(([0,0,0],))
51    elif SGData['SGLatt'] == 'A':
52        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5]))
53    elif SGData['SGLatt'] == 'B':
54        SGData['SGCen'] = np.array(([0,0,0],[.5,0,.5]))
55    elif SGData['SGLatt'] == 'C':
56        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,0,]))
57    elif SGData['SGLatt'] == 'I':
58        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,.5]))
59    elif SGData['SGLatt'] == 'F':
60        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5],[.5,0,.5],[.5,.5,0,]))
61    elif SGData['SGLatt'] == 'R':
62        SGData['SGCen'] = np.array(([0,0,0],[1./3.,2./3.,2./3.],[2./3.,1./3.,1./3.]))
63    SGData['SGOps'] = []
64    for i in range(SGInfo[5]):
65        Mat = np.array(SGInfo[6][i])
66        Trns = np.array(SGInfo[7][i])
67        SGData['SGOps'].append([Mat,Trns])
68    if SGData['SGLaue'] in '-1':
69        SGData['SGSys'] = SysSym[0]
70    elif SGData['SGLaue'] in '2/m':
71        SGData['SGSys'] = SysSym[1]
72    elif SGData['SGLaue'] in 'mmm':
73        SGData['SGSys'] = SysSym[2]
74    elif SGData['SGLaue'] in ['4/m','4/mmm']:
75        SGData['SGSys'] = SysSym[3]
76    elif SGData['SGLaue'] in ['3R','3mR']:
77        SGData['SGSys'] = SysSym[4]
78    elif SGData['SGLaue'] in ['3','3m1','31m']:
79        SGData['SGSys'] = SysSym[5]
80    elif SGData['SGLaue'] in ['6/m','6/mmm']:
81        SGData['SGSys'] = SysSym[6]
82    elif SGData['SGLaue'] in ['m3','m3m']:
83        SGData['SGSys'] = SysSym[7]
84    SGData['SGPolax'] = SGpolar(SGData)
85    return SGInfo[8],SGData
86
87def SGErrors(IErr):
88    '''
89    Interprets the error message code from SpcGroup. Used in SpaceGroup.
90    input:
91        SGError - from SpcGroup
92    returns:
93        ErrString - a string with the error message or "Unknown error"
94    '''
95
96    ErrString = [' ',
97        'Less than 2 operator fields were found',
98        'Illegal Lattice type, not P, A, B, C, I, F or R',
99        'Rhombohedral lattice requires a 3-axis',
100        'Minus sign does not preceed 1, 2, 3, 4 or 6',
101        'Either a 5-axis anywhere or a 3-axis in field not allowed',
102        ' ',
103        'I for COMPUTED GO TO out of range.',
104        'An a-glide mirror normal to A not allowed',
105        'A b-glide mirror normal to B not allowed',
106        'A c-glide mirror normal to C not allowed',
107        'D-glide in a primitive lattice not allowed',
108        'A 4-axis not allowed in the 2nd operator field',
109        'A 6-axis not allowed in the 2nd operator field',
110        'More than 24 matrices needed to define group',
111        ' ',
112        'Improper construction of a rotation operator',
113        'Mirror following a / not allowed',
114        'A translation conflict between operators',
115        'The 2bar operator is not allowed',
116        '3 fields are legal only in R & m3 cubic groups',
117        'Syntax error. Expected I -4 3 d at this point',
118        ' ',
119        'A or B centered tetragonal not allowed',
120        ' ','unknown error in sgroup',' ',' ',' ',
121        'Illegal character in the space group symbol',
122        ]
123    try:
124        return ErrString[IErr]
125    except:
126        return "Unknown error"
127def SGpolar(SGData):
128    '''
129    Determine identity of polar axes if any
130    '''
131    POL = ('','x','y','x y','z','x z','y z','xyz','111')
132    NP = [1,2,4]
133    NPZ = [0,1]
134    for M,T in SGData['SGOps']:
135        for i in range(3):
136            if M[i][i] <= 0.: NP[i] = 0
137        if M[0][2] > 0: NPZ[0] = 8
138        if M[1][2] > 0: NPZ[1] = 0
139    NPol = (NP[0]+NP[1]+NP[2]+NPZ[0]*NPZ[1])*(1-int(SGData['SGInv']))
140    return POL[NPol]
141   
142def SGPrint(SGData):
143    '''
144    Print the output of SpcGroup in a nicely formatted way. Used in SpaceGroup
145    input:
146        SGData - from SpcGroup
147    returns:
148        SGText - list of strings with the space group details
149    '''
150    Mult = len(SGData['SGCen'])*len(SGData['SGOps'])*(int(SGData['SGInv'])+1)
151    SGText = []
152    SGText.append(' Space Group: '+SGData['SpGrp'])
153    CentStr = 'centrosymmetric'
154    if not SGData['SGInv']:
155        CentStr = 'non'+CentStr
156    if SGData['SGLatt'] in 'ABCIFR':
157        SGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
158    else:
159        SGText.append(' The lattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())       
160    SGText.append(' Multiplicity of a general site is '+str(Mult))
161    SGText.append(' The Laue symmetry is '+SGData['SGLaue'])
162    if SGData['SGUniq'] in ['a','b','c']:
163        SGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
164    if SGData['SGInv']:
165        SGText.append(' The inversion center is located at 0,0,0')
166    if SGData['SGPolax']:
167        SGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
168    SGText.append('\n'+' The equivalent positions are:')
169    if SGData['SGLatt'] != 'P':
170        SGText.append('\n ('+Latt2text(SGData['SGLatt'])+')+\n')
171    Ncol = 2
172    line = ' '
173    col = 0
174    for iop,[M,T] in enumerate(SGData['SGOps']):
175        OPtxt = MT2text(M,T)
176        Fld = '(%2i) '%(iop+1)+OPtxt+'\t'
177        line += Fld
178        if '/' not in Fld:
179            line += '\t'
180        col += 1
181        if col == Ncol:
182            SGText.append(line)       
183            line = ' '
184            col = 0
185    SGText.append(line)       
186    return SGText
187   
188def MT2text(M,T):
189    #From space group matrix/translation operator returns text version
190    XYZ = ('-Z ','-Y ','-X ','X-Y','ERR','Y-X',' X ',' Y ',' Z ','+X ','+Y ','+Z ')
191    TRA = ('   ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
192    Fld = ''
193    for j in range(3):
194        IJ = int(round(2*M[j][0]+3*M[j][1]+4*M[j][2]+4))%12
195        IK = int(round(T[j]*12))%12
196        if IK > 0 and IJ > 4: IJ += 3
197        Fld += TRA[IK]+XYZ[IJ]
198        if j != 2: Fld += ','
199    return Fld
200   
201def Latt2text(Latt):
202    #From lattice type ('P',A', etc.) returns ';' delimited cell centering vectors
203    lattTxt = {'A':'0,0,0; 0,1/2,1/2','B':'0,0,0; 1/2,0,1/2',
204        'C':'0,0,0; 1/2,1/2,0','I':'0,0,0; 1/2,1/2,1/2',
205        'F':'0,0,0; 0,1/2,1/2; 1/2,0,1/2; 1/2,1/2,0',
206        'R':'0,0,0; 1/3,2/3,2/3; 2/3,1/3,1/3','P':'0,0,0'}
207    return lattTxt[Latt]   
208       
209def SpaceGroup(SGSymbol):
210    '''
211    Print the output of SpcGroup in a nicely formatted way.
212    input:
213        SGSymbol - space group symbol (string) with spaces between axial fields
214    returns:
215        nothing
216    '''
217    E,A = SpcGroup(SGSymbol)
218    if E > 0:
219        print SGErrors(E)
220        return
221    for l in SGPrint(A):
222        print l
223
224def MoveToUnitCell(xyz):
225    '''
226    Translates a set of coordinates so that all values are >=0 and < 1
227    input:
228        xyz - a list or numpy array of fractional coordinates
229    returns:
230        XYZ - numpy array of new coordinates inside 0-1
231    '''
232    XYZ = np.zeros(3)
233    for i,x in enumerate(xyz):
234        XYZ[i] = (x-int(x))%1.0
235    return XYZ
236       
237def Opposite(XYZ,toler=0.0002):
238    '''
239    Gives opposite corner, edge or face of unit cell for position within tolerance.
240        Result may be just outside the cell within tolerance
241    input:
242        XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
243        toler: unit cell fraction tolerance making opposite
244    returns:
245        XYZ: array of opposite positions; always contains XYZ
246    '''
247    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]]
248    TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
249    perm = TB*perm3
250    cperm = ['%d%d%d'%(i,j,k) for i,j,k in perm]
251    D = dict(zip(cperm,perm))
252    new = []
253    for key in D:
254        new.append(np.array(D[key])+np.array(XYZ))
255    return new
256       
257def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
258    '''
259    Generates the equivalent positions for a specified coordinate and space group
260    input: 
261        XYZ an array, tuple or list containing 3 elements: x, y & z
262        SGData, from SpcGroup
263        All  = True return all equivalent positions including duplicates
264             = False return only unique positions
265        Uij  = [U11,U22,U33,U12,U13,U23] or [] if no Uij
266        Move = True move generated atom positions to be inside cell
267             = False do not move atoms       
268    return: [[XYZEquiv],Idup,[UijEquiv]]
269        [XYZEquiv] is list of equivalent positions (XYZ is first entry)
270        Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
271        is centering operator number (1-4) and - is for inversion
272        Cell = unit cell translations needed to put new positions inside cell
273        [UijEquiv] - equivalent Uij; absent if no Uij given
274    '''
275    XYZEquiv = []
276    UijEquiv = []
277    Idup = []
278    Cell = []
279    X = np.array(XYZ)
280    if Move:
281        X = MoveToUnitCell(X)
282    for ic,cen in enumerate(SGData['SGCen']):
283        C = np.array(cen)
284        for invers in range(int(SGData['SGInv']+1)):
285            for io,[M,T] in enumerate(SGData['SGOps']):
286                idup = ((io+1)+100*ic)*(1-2*invers)
287                XT = np.inner(M,X)+T
288                if len(Uij):
289                    U = Uij2U(Uij)
290                    U = np.inner(M,np.inner(U,M).T)
291                    newUij = U2Uij(U)
292                if invers:
293                    XT = -XT
294                XT += C
295                newX = MoveToUnitCell(XT)
296                cell = np.asarray(np.rint(newX-XT),dtype=np.int32)
297                if All:
298                    if np.allclose(newX,X,atol=0.0002):
299                        idup = False
300                else:
301                    if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
302                        idup = False
303                if All or idup:
304                    XYZEquiv.append(newX)
305                    Idup.append(idup)
306                    Cell.append(cell)
307                    if len(Uij):
308                        UijEquiv.append(newUij)                   
309    if len(Uij):
310        return zip(XYZEquiv,UijEquiv,Idup,Cell)
311    else:
312        return zip(XYZEquiv,Idup,Cell)
313
314def GenHKLf(HKL,SGData,Friedel=False):
315    '''
316    Uses old GSAS Fortran routine genhkl.for
317    input:
318        HKL - [h,k,l]
319        SGData - space group data obtained from SpcGroup
320        Friedel = True to retain Friedel pairs for centrosymmetric case
321    returns:
322        iabsnt = True is reflection is forbidden by symmetry
323        mulp = reflection multiplicity including Fridel pairs
324        Uniq = numpy array of equivalent hkl in descending order of h,k,l
325    '''
326    hklf = HKL+[0,]
327    Ops = SGData['SGOps']
328    OpM = np.array([op[0] for op in Ops])
329    OpT = np.array([op[1] for op in Ops])
330    Inv = SGData['SGInv']
331    Cen = np.array([cen for cen in SGData['SGCen']])
332   
333    Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
334    h,k,l,f = Uniq
335    Uniq=np.array(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq]))
336    phi = f[:Nuniq]
337    Uniq = np.array(Uniq)
338   
339    return iabsnt,2*mulp,Uniq,phi       #include Friedel pairs in powder mulp
340                                 
341def GetOprPtrName(key):           
342    OprPtrName = {
343        '-6643':[   2,' 1bar ', 1],'6479' :[  10,'  2z  ', 2],'-6479':[   9,'  mz  ', 3],
344        '6481' :[   7,'  my  ', 4],'-6481':[   6,'  2y  ', 5],'6641' :[   4,'  mx  ', 6],
345        '-6641':[   3,'  2x  ', 7],'6591' :[  28,' m+-0 ', 8],'-6591':[  27,' 2+-0 ', 9],
346        '6531' :[  25,' m110 ',10],'-6531':[  24,' 2110 ',11],'6537' :[  61,'  4z  ',12],
347        '-6537':[  62,' -4z  ',13],'975'  :[  68,' 3+++1',14],'6456' :[ 114,'  3z1 ',15],
348        '-489' :[  73,' 3+-- ',16],'483'  :[  78,' 3-+- ',17],'-969' :[  83,' 3--+ ',18],
349        '819'  :[  22,' m+0- ',19],'-819' :[  21,' 2+0- ',20],'2431' :[  16,' m0+- ',21],
350        '-2431':[  15,' 20+- ',22],'-657' :[  19,' m101 ',23],'657'  :[  18,' 2101 ',24],
351        '1943' :[  48,' -4x  ',25],'-1943':[  47,'  4x  ',26],'-2429':[  13,' m011 ',27],
352        '2429' :[  12,' 2011 ',28],'639'  :[  55,' -4y  ',29],'-639' :[  54,'  4y  ',30],
353        '-6484':[ 146,' 2010 ', 4],'6484' :[ 139,' m010 ', 5],'-6668':[ 145,' 2100 ', 6],
354        '6668' :[ 138,' m100 ', 7],'-6454':[ 148,' 2120 ',18],'6454' :[ 141,' m120 ',19],
355        '-6638':[ 149,' 2210 ',20],'6638' :[ 142,' m210 ',21],              #search ends here
356        '2223' :[  68,' 3+++2',39],
357        '6538' :[ 106,'  6z1 ',40],'-2169':[  83,' 3--+2',41],'2151' :[  73,' 3+--2',42],
358        '2205' :[  79,'-3-+-2',43],'-2205':[  78,' 3-+-2',44],'489'  :[  74,'-3+--1',45],
359        '801'  :[  53,'  4y1 ',46],'1945' :[  47,'  4x3 ',47],'-6585':[  62,' -4z3 ',48],
360        '6585' :[  61,'  4z3 ',49],'6584' :[ 114,'  3z2 ',50],'6666' :[ 106,'  6z5 ',51],
361        '6643' :[   1,' Iden ',52],'-801' :[  55,' -4y1 ',53],'-1945':[  48,' -4x3 ',54],
362        '-6666':[ 105,' -6z5 ',55],'-6538':[ 105,' -6z1 ',56],'-2223':[  69,'-3+++2',57],
363        '-975' :[  69,'-3+++1',58],'-6456':[ 113,' -3z1 ',59],'-483' :[  79,'-3-+-1',60],
364        '969'  :[  84,'-3--+1',61],'-6584':[ 113,' -3z2 ',62],'2169' :[  84,'-3--+2',63],
365        '-2151':[  74,'-3+--2',64],'0':[0,' ????',0]
366        }
367    return OprPtrName[key]
368
369def GetKNsym(key):
370    KNsym = {
371        '0'         :'    1   ','1'         :'   -1   ','64'        :'  2(100)','32'        :'  m(100)',
372        '97'        :'2/m(100)','16'        :'  2(010)','8'         :'  m(010)','25'        :'2/m(010)',
373        '2'         :'  2(001)','4'         :'  m(001)','7'         :'2/m(001)','134217728' :'  2(011)',
374        '67108864'  :'  m(011)','201326593' :'2/m(011)','2097152'   :'  2(0+-)','1048576'   :'  m(0+-)',
375        '3145729'   :'2/m(0+-)','8388608'   :'  2(101)','4194304'   :'  m(101)','12582913'  :'2/m(101)',
376        '524288'    :'  2(+0-)','262144'    :'  m(+0-)','796433'    :'2/m(+0-)','1024'      :'  2(110)',
377        '512'       :'  m(110)','1537'      :'2/m(110)','256'       :'  2(+-0)','128'       :'  m(+-0)',
378        '385'       :'2/m(+-0)','76'        :'mm2(100)','52'        :'mm2(010)','42'        :'mm2(001)',
379        '135266336' :'mm2(011)','69206048'  :'mm2(0+-)','8650760'   :'mm2(101)','4718600'   :'mm2(+0-)',
380        '1156'      :'mm2(110)','772'       :'mm2(+-0)','82'        :'  222   ','136314944' :'222(100)',
381        '8912912'   :'222(010)','1282'      :'222(001)','127'       :'  mmm   ','204472417' :'mmm(100)',
382        '13369369'  :'mmm(010)','1927'      :'mmm(001)','33554496'  :'  4(100)','16777280'  :' -4(100)',
383        '50331745'  :'4/m(100)','169869394' :'422(100)','84934738'  :'-42m 100','101711948' :'4mm(100)',
384        '254804095' :'4/mmm100','536870928 ':'  4(010)','268435472' :' -4(010)','805306393' :'4/m (10)',
385        '545783890' :'422(010)','272891986' :'-42m 010','541327412' :'4mm(010)','818675839' :'4/mmm010',
386        '2050'      :'  4(001)','4098'      :' -4(001)','6151'      :'4/m(001)','3410'      :'422(001)',
387        '4818'      :'-42m 001','2730'      :'4mm(001)','8191'      :'4/mmm001','8192'      :'  3(111)',
388        '8193'      :' -3(111)','2629888'   :' 32(111)','1319040'   :' 3m(111)','3940737'   :'-3m(111)',
389        '32768'     :'  3(+--)','32769'     :' -3(+--)','10519552'  :' 32(+--)','5276160'   :' 3m(+--)',
390        '15762945'  :'-3m(+--)','65536'     :'  3(-+-)','65537'     :' -3(-+-)','134808576' :' 32(-+-)',
391        '67437056'  :' 3m(-+-)','202180097' :'-3m(-+-)','131072'    :'  3(--+)','131073'    :' -3(--+)',
392        '142737664' :' 32(--+)','71434368'  :' 3m(--+)','214040961' :'-3m(--+)','237650'    :'   23   ',
393        '237695'    :'   m3   ','715894098' :'   432  ','358068946' :'  -43m  ','1073725439':'   m3m  ',
394        '68157504'  :' mm2d100','4456464'   :' mm2d010','642'       :' mm2d001','153092172' :'-4m2 100',
395        '277348404' :'-4m2 010','5418'      :'-4m2 001','1075726335':'  6/mmm ','1074414420':'-6m2 100',
396        '1075070124':'-6m2 120','1075069650':'   6mm  ','1074414890':'   622  ','1073758215':'   6/m  ',
397        '1073758212':'   -6   ','1073758210':'    6   ','1073759865':'-3m(100)','1075724673':'-3m(120)',
398        '1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
399        '1073758209':'   -3   ','1073758208':'    3   ','1074135143':'mmm(100)','1075314719':'mmm(010)',
400        '1073743751':'mmm(110)','1074004034':' mm2z100','1074790418':' mm2z010','1073742466':' mm2z110',
401        '1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
402        '1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
403        '1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
404        '1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
405        '1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
406        '1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
407        '1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
408        '1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':'   -1   '
409        }
410    return KNsym[key]       
411
412def GetNXUPQsym(siteSym):       
413    NXUPQsym = {
414        '    1   ':(28,29,28,28),'   -1   ':( 1,29,28, 0),'  2(100)':(12,18,12,25),'  m(100)':(25,18,12,25),
415        '2/m(100)':( 1,18, 0,-1),'  2(010)':(13,17,13,24),'  m(010)':(24,17,13,24),'2/m(010)':( 1,17, 0,-1),
416        '  2(001)':(14,16,14,23),'  m(001)':(23,16,14,23),'2/m(001)':( 1,16, 0,-1),'  2(011)':(10,23,10,22),
417        '  m(011)':(22,23,10,22),'2/m(011)':( 1,23, 0,-1),'  2(0+-)':(11,24,11,21),'  m(0+-)':(21,24,11,21),
418        '2/m(0+-)':( 1,24, 0,-1),'  2(101)':( 8,21, 8,20),'  m(101)':(20,21, 8,20),'2/m(101)':( 1,21, 0,-1),
419        '  2(+0-)':( 9,22, 9,19),'  m(+0-)':(19,22, 9,19),'2/m(+0-)':( 1,22, 0,-1),'  2(110)':( 6,19, 6,18),
420        '  m(110)':(18,19, 6,18),'2/m(110)':( 1,19, 0,-1),'  2(+-0)':( 7,20, 7,17),'  m(+-0)':(17,20, 7,17),
421        '2/m(+-0)':( 1,20, 0,-1),'mm2(100)':(12,10, 0,-1),'mm2(010)':(13,10, 0,-1),'mm2(001)':(14,10, 0,-1),
422        'mm2(011)':(10,13, 0,-1),'mm2(0+-)':(11,13, 0,-1),'mm2(101)':( 8,12, 0,-1),'mm2(+0-)':( 9,12, 0,-1),
423        'mm2(110)':( 6,11, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'  222   ':( 1,10, 0,-1),'222(100)':( 1,13, 0,-1),
424        '222(010)':( 1,12, 0,-1),'222(001)':( 1,11, 0,-1),'  mmm   ':( 1,10, 0,-1),'mmm(100)':( 1,13, 0,-1),
425        'mmm(010)':( 1,12, 0,-1),'mmm(001)':( 1,11, 0,-1),'  4(100)':(12, 4,12, 0),' -4(100)':( 1, 4,12, 0),
426        '4/m(100)':( 1, 4,12,-1),'422(100)':( 1, 4, 0,-1),'-42m 100':( 1, 4, 0,-1),'4mm(100)':(12, 4, 0,-1),
427        '4/mmm100':( 1, 4, 0,-1),'  4(010)':(13, 3,13, 0),' -4(010)':( 1, 3,13, 0),'4/m (10)':( 1, 3,13,-1),
428        '422(010)':( 1, 3, 0,-1),'-42m 010':( 1, 3, 0,-1),'4mm(010)':(13, 3, 0,-1),'4/mmm010':(1, 3, 0,-1,),
429        '  4(001)':(14, 2,14, 0),' -4(001)':( 1, 2,14, 0),'4/m(001)':( 1, 2,14,-1),'422(001)':( 1, 2, 0,-1),
430        '-42m 001':( 1, 2, 0,-1),'4mm(001)':(14, 2, 0,-1),'4/mmm001':( 1, 2, 0,-1),'  3(111)':( 2, 5, 2, 0),
431        ' -3(111)':( 1, 5, 2, 0),' 32(111)':( 1, 5, 0, 2),' 3m(111)':( 2, 5, 0, 2),'-3m(111)':( 1, 5, 0,-1),
432        '  3(+--)':( 5, 8, 5, 0),' -3(+--)':( 1, 8, 5, 0),' 32(+--)':( 1, 8, 0, 5),' 3m(+--)':( 5, 8, 0, 5),
433        '-3m(+--)':( 1, 8, 0,-1),'  3(-+-)':( 4, 7, 4, 0),' -3(-+-)':( 1, 7, 4, 0),' 32(-+-)':( 1, 7, 0, 4),
434        ' 3m(-+-)':( 4, 7, 0, 4),'-3m(-+-)':( 1, 7, 0,-1),'  3(--+)':( 3, 6, 3, 0),' -3(--+)':( 1, 6, 3, 0),
435        ' 32(--+)':( 1, 6, 0, 3),' 3m(--+)':( 3, 6, 0, 3),'-3m(--+)':( 1, 6, 0,-1),'   23   ':( 1, 1, 0, 0),
436        '   m3   ':( 1, 1, 0, 0),'   432  ':( 1, 1, 0, 0),'  -43m  ':( 1, 1, 0, 0),'   m3m  ':( 1, 1, 0, 0),
437        ' mm2d100':(12,13, 0,-1),' mm2d010':(13,12, 0,-1),' mm2d001':(14,11, 0,-1),'-4m2 100':( 1, 4, 0,-1),
438        '-4m2 010':( 1, 3, 0,-1),'-4m2 001':( 1, 2, 0,-1),'  6/mmm ':( 1, 9, 0,-1),'-6m2 100':( 1, 9, 0,-1),
439        '-6m2 120':( 1, 9, 0,-1),'   6mm  ':(14, 9, 0,-1),'   622  ':( 1, 9, 0,-1),'   6/m  ':( 1, 9,14,-1),
440        '   -6   ':( 1, 9,14, 0),'    6   ':(14, 9,14, 0),'-3m(100)':( 1, 9, 0,-1),'-3m(120)':( 1, 9, 0,-1),
441        ' 3m(100)':(14, 9, 0,14),' 3m(120)':(14, 9, 0,14),' 32(100)':( 1, 9, 0,14),' 32(120)':( 1, 9, 0,14),
442        '   -3   ':( 1, 9,14, 0),'    3   ':(14, 9,14, 0),'mmm(100)':( 1,14, 0,-1),'mmm(010)':( 1,15, 0,-1),
443        'mmm(110)':( 1,11, 0,-1),' mm2z100':(14,14, 0,-1),' mm2z010':(14,15, 0,-1),' mm2z110':(14,11, 0,-1),
444        'mm2(100)':(12,14, 0,-1),'mm2(010)':(13,15, 0,-1),'mm2(110)':( 6,11, 0,-1),'mm2(120)':(15,14, 0,-1),
445        'mm2(210)':(16,15, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'222(100)':( 1,14, 0,-1),'222(010)':( 1,15, 0,-1),
446        '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),
447        '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),
448        ' m(001) ':(23,16,14,23),' m(100) ':(26,25,12,26),' m(010) ':(27,28,13,27),' m(110) ':(18,19, 6,18),
449        ' m(120) ':(24,27,15,24),' m(210) ':(25,26,16,25),' m(+-0) ':(17,20, 7,17),' 2(001) ':(14,16,14,23),
450        ' 2(100) ':(12,25,12,26),' 2(010) ':(13,28,13,27),' 2(110) ':( 6,19, 6,18),' 2(120) ':(15,27,15,24),
451        ' 2(210) ':(16,26,16,25),' 2(+-0) ':( 7,20, 7,17),'   -1   ':( 1,29,28, 0)
452        }
453    return NXUPQsym[siteSym]
454
455def GetCSxinel(siteSym): 
456    CSxinel = [[],                         # 0th empty - indices are Fortran style
457        [[0,0,0],[ 0.0, 0.0, 0.0]],      #  0  0  0
458        [[1,1,1],[ 1.0, 1.0, 1.0]],      #  X  X  X
459        [[1,1,1],[ 1.0, 1.0,-1.0]],      #  X  X -X
460        [[1,1,1],[ 1.0,-1.0, 1.0]],      #  X -X  X
461        [[1,1,1],[ 1.0,-1.0,-1.0]],      # -X  X  X
462        [[1,1,0],[ 1.0, 1.0, 0.0]],      #  X  X  0
463        [[1,1,0],[ 1.0,-1.0, 0.0]],      #  X -X  0
464        [[1,0,1],[ 1.0, 0.0, 1.0]],      #  X  0  X
465        [[1,0,1],[ 1.0, 0.0,-1.0]],      #  X  0 -X
466        [[0,1,1],[ 0.0, 1.0, 1.0]],      #  0  Y  Y
467        [[0,1,1],[ 0.0, 1.0,-1.0]],      #  0  Y -Y
468        [[1,0,0],[ 1.0, 0.0, 0.0]],      #  X  0  0
469        [[0,1,0],[ 0.0, 1.0, 0.0]],      #  0  Y  0
470        [[0,0,1],[ 0.0, 0.0, 1.0]],      #  0  0  Z
471        [[1,1,0],[ 1.0, 2.0, 0.0]],      #  X 2X  0
472        [[1,1,0],[ 2.0, 1.0, 0.0]],      # 2X  X  0
473        [[1,1,2],[ 1.0, 1.0, 1.0]],      #  X  X  Z
474        [[1,1,2],[ 1.0,-1.0, 1.0]],      #  X -X  Z
475        [[1,2,1],[ 1.0, 1.0, 1.0]],      #  X  Y  X
476        [[1,2,1],[ 1.0, 1.0,-1.0]],      #  X  Y -X
477        [[1,2,2],[ 1.0, 1.0, 1.0]],      #  X  Y  Y
478        [[1,2,2],[ 1.0, 1.0,-1.0]],      #  X  Y -Y
479        [[1,2,0],[ 1.0, 1.0, 0.0]],      #  X  Y  0
480        [[1,0,2],[ 1.0, 0.0, 1.0]],      #  X  0  Z
481        [[0,1,2],[ 0.0, 1.0, 1.0]],      #  0  Y  Z
482        [[1,1,2],[ 1.0, 2.0, 1.0]],      #  X 2X  Z
483        [[1,1,2],[ 2.0, 1.0, 1.0]],      # 2X  X  Z
484        [[1,2,3],[ 1.0, 1.0, 1.0]],      #  X  Y  Z
485        ]
486    indx = GetNXUPQsym(siteSym)
487    return CSxinel[indx[0]]
488   
489def GetCSuinel(siteSym):
490    # returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers
491    CSuinel = [[],                                             # 0th empty - indices are Fortran style
492        [[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]],    #  A  A  A  0  0  0
493        [[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]],    #  A  A  C  0  0  0
494        [[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]],    #  A  B  A  0  0  0
495        [[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]],    #  A  B  B  0  0  0
496        [[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]],    #  A  A  A  D  D  D
497        [[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]],    #  A  A  A  D -D -D
498        [[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]],    #  A  A  A  D -D  D
499        [[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]],    #  A  A  A  D  D -D
500        [[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]],    #  A  A  C A/2 0  0
501        [[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]],    #  A  B  C  0  0  0
502        [[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]],    #  A  A  C  D  0  0
503        [[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]],    #  A  B  A  0  E  0
504        [[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]],    #  A  B  B  0  0  F
505        [[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]],    #  A  B  C B/2 0  0
506        [[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]],    #  A  B  C A/2 0  0
507        [[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]],    #  A  B  C  D  0  0
508        [[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]],    #  A  B  C  0  E  0
509        [[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]],    #  A  B  C  0  0  F
510        [[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]],    #  A  A  C  D  E -E
511        [[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]],    #  A  A  C  D  E  E
512        [[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]],    #  A  B  A  D  E -D
513        [[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]],    #  A  B  A  D  E  D
514        [[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]],    #  A  B  B  D -D  F
515        [[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]],    #  A  B  B  D  D  F
516        [[1,2,3,2,4,4],[ 1.0, 1.0, 1.0, 0.5, 0.5, 1.0],[1,1,1,0,0,1],[1.0,1.0,1.0,0.5,0.0,0.0]],    #  A  B  C B/2 F/2 F
517        [[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]],    #  A  B  C A/2  0  F
518        [[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]],    #  A  B  C B/2  E  0
519        [[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]],    #  A  B  C A/2  E E/2
520        [[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]],    #  A  B  C  D  E   F
521        ]
522    indx = GetNXUPQsym(siteSym)
523    return CSuinel[indx[1]]
524   
525def MustrainNames(SGData):
526    laue = SGData['SGLaue']
527    uniq = SGData['SGUniq']
528    if laue in ['m3','m3m']:
529        return ['S400','S220']
530    elif laue in ['6/m','6/mmm','3m1']:
531        return ['S400','S004','S202']
532    elif laue in ['31m','3']:
533        return ['S400','S004','S202','S211']
534    elif laue in ['3R','3mR']:
535        return ['S400','S220','S310','S211']
536    elif laue in ['4/m','4/mmm']:
537        return ['S400','S004','S220','S022']
538    elif laue in ['mmm']:
539        return ['S400','S040','S004','S220','S202','S022']
540    elif laue in ['2/m']:
541        SHKL = ['S400','S040','S004','S220','S202','S022']
542        if uniq == 'a':
543            SHKL += ['S013','S031','S211']
544        elif uniq == 'b':
545            SHKL += ['S301','S103','S121']
546        elif uniq == 'c':
547            SHKL += ['S130','S310','S112']
548        return SHKL
549    else:
550        SHKL = ['S400','S040','S004','S220','S202','S022']
551        SHKL += ['S310','S103','S031','S130','S301','S013']
552        SHKL += ['S211','S121','S112']
553        return SHKL
554
555def HStrainNames(SGData):
556    laue = SGData['SGLaue']
557    uniq = SGData['SGUniq']
558    if laue in ['m3','m3m']:
559        return ['D11','eA']         #add cubic strain term
560    elif laue in ['6/m','6/mmm','3m1','31m','3']:
561        return ['D11','D33']
562    elif laue in ['3R','3mR']:
563        return ['D11','D12']
564    elif laue in ['4/m','4/mmm']:
565        return ['D11','D33']
566    elif laue in ['mmm']:
567        return ['D11','D22','D33']
568    elif laue in ['2/m']:
569        Dij = ['D11','D22','D33']
570        if uniq == 'a':
571            Dij += ['D23']
572        elif uniq == 'b':
573            Dij += ['D13']
574        elif uniq == 'c':
575            Dij += ['D12']
576        return Dij
577    else:
578        Dij = ['D11','D22','D33','D12','D13','D23']
579        return Dij
580   
581def MustrainCoeff(HKL,SGData):
582    #NB: order of terms is the same as returned by MustrainNames
583    laue = SGData['SGLaue']
584    uniq = SGData['SGUniq']
585    h,k,l = HKL
586    Strm = []
587    if laue in ['m3','m3m']:
588        Strm.append(h**4+k**4+l**4)
589        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
590    elif laue in ['6/m','6/mmm','3m1']:
591        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
592        Strm.append(l**4)
593        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
594    elif laue in ['31m','3']:
595        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
596        Strm.append(l**4)
597        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
598        Strm.append(4.0*h*k*l*(h+k))
599    elif laue in ['3R','3mR']:
600        Strm.append(h**4+k**4+l**4)
601        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
602        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
603        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
604    elif laue in ['4/m','4/mmm']:
605        Strm.append(h**4+k**4)
606        Strm.append(l**4)
607        Strm.append(3.0*(h*k)**2)
608        Strm.append(3.0*((h*l)**2+(k*l)**2))
609    elif laue in ['mmm']:
610        Strm.append(h**4)
611        Strm.append(k**4)
612        Strm.append(l**4)
613        Strm.append(3.0*(h*k)**2)
614        Strm.append(3.0*(h*l)**2)
615        Strm.append(3.0*(k*l)**2)
616    elif laue in ['2/m']:
617        Strm.append(h**4)
618        Strm.append(k**4)
619        Strm.append(l**4)
620        Strm.append(3.0*(h*k)**2)
621        Strm.append(3.0*(h*l)**2)
622        Strm.append(3.0*(k*l)**2)
623        if uniq == 'a':
624            Strm.append(2.0*k*l**3)
625            Strm.append(2.0*l*k**3)
626            Strm.append(4.0*k*l*h**2)
627        elif uniq == 'b':
628            Strm.append(2.0*l*h**3)
629            Strm.append(2.0*h*l**3)
630            Strm.append(4.0*h*l*k**2)
631        elif uniq == 'c':
632            Strm.append(2.0*h*k**3)
633            Strm.append(2.0*k*h**3)
634            Strm.append(4.0*h*k*l**2)
635    else:
636        Strm.append(h**4)
637        Strm.append(k**4)
638        Strm.append(l**4)
639        Strm.append(3.0*(h*k)**2)
640        Strm.append(3.0*(h*l)**2)
641        Strm.append(3.0*(k*l)**2)
642        Strm.append(2.0*k*h**3)
643        Strm.append(2.0*h*l**3)
644        Strm.append(2.0*l*k**3)
645        Strm.append(2.0*h*k**3)
646        Strm.append(2.0*l*h**3)
647        Strm.append(2.0*k*l**3)
648        Strm.append(4.0*k*l*h**2)
649        Strm.append(4.0*h*l*k**2)
650        Strm.append(4.0*k*h*l**2)
651    return Strm
652   
653def Muiso2Shkl(muiso,SGData,cell):
654    #this is to convert isotropic mustrain to generalized Shkls - doesn't work just now
655    import GSASIIlattice as G2lat
656    from scipy.optimize import fmin
657    A = G2lat.cell2AB(cell)[0]
658    def minMus(Shkl,H,muiso,SGData,A):
659        U = np.inner(A.T,H)
660        S = np.array(MustrainCoeff(H.T,SGData))
661        sum = np.sqrt(np.sum(np.multiply(S,Shkl)))
662        return abs(muiso-sum*H)
663    laue = SGData['SGLaue']
664    if laue in ['m3','m3m']:
665        H = [[1,0,0],[1,1,0]]
666        S0 = [0.01,0.01]
667    elif laue in ['6/m','6/mmm','3m1']:
668        H = [[1,0,0],[0,0,1],[1,0,1]]
669        S0 = [0.01,0.01,0.01]
670    elif laue in ['31m','3']:
671        H = [[1,0,0],[0,0,1],[1,0,1],[1,1,1]]
672        S0 = [0.01,0.01,0.01,0.01]
673    elif laue in ['3R','3mR']:
674        H = [[1,0,0],[1,1,0],[1,0,1],[1,1,1]]
675        S0 = [0.01,0.01,0.01,0.01]
676    elif laue in ['4/m','4/mmm']:
677        H = [[1,0,0],[0,0,1],[1,1,0],[1,0,1]]
678        S0 = [0.01,0.01,0.01,0.01]
679    elif laue in ['mmm']:
680        H = [[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,0,1],[0,1,1]]
681        S0 = [0.01,0.01,0.01,0.01,0.01,0.01]
682    elif laue in ['2/m']:
683        H = [[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,0,1],[0,1,1]]
684        if uniq == 'a':
685            H.append([0,1,-1])
686            H.append([0,-2,1])
687        elif uniq == 'b':
688            H.append([1,0,-1])
689            H.append([-2,0,1])
690        elif uniq == 'c':
691            H.append([1,-1,0])
692            H.append([-2,1,0])
693        H.append([1,1,1])
694        S0 = [9*[0.01,]]
695    else:
696        H = [[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,0,1],[0,1,1],
697            [-1,1,0],[1,0,-1],[0,-1,1],[1,-2,0],[-2,0,1],[0,1,-2],
698            [1,-1,1],[-1, 1, 1],[1,-1,1]]
699        S0 = [15*[0.01,]]
700    H = np.array(H)
701    S0 = np.array(S0)
702    return fmin(minMus,S0,(H,muiso,SGData,A))
703       
704def SytSym(XYZ,SGData):
705    '''
706    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
707    input: 
708       XYZ: an array, tuple or list containing 3 elements: x, y & z
709       SGData: from SpcGroup
710    Returns a two element tuple:
711       The 1st element is a code for the site symmetry (see GetKNsym)
712       The 2nd element is the site multiplicity
713    '''
714    def PackRot(SGOps):
715        IRT = []
716        for ops in SGOps:
717            M = ops[0]
718            irt = 0
719            for j in range(2,-1,-1):
720                for k in range(2,-1,-1):
721                    irt *= 3
722                    irt += M[k][j]
723            IRT.append(int(irt))
724        return IRT
725       
726    SymName = ''
727    Mult = 1
728    Isym = 0
729    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
730        Isym = 1073741824
731    Jdup = 0
732    Xeqv = GenAtom(XYZ,SGData,True)
733    IRT = PackRot(SGData['SGOps'])
734    L = -1
735    for ic,cen in enumerate(SGData['SGCen']):
736        for invers in range(int(SGData['SGInv']+1)):
737            for io,ops in enumerate(SGData['SGOps']):
738                irtx = (1-2*invers)*IRT[io]
739                L += 1
740                if not Xeqv[L][1]:
741                    Jdup += 1
742                    jx = GetOprPtrName(str(irtx))
743                    if jx[2] < 39:
744                        Isym += 2**(jx[2]-1)
745    if Isym == 1073741824: Isym = 0
746    Mult = len(SGData['SGOps'])*len(SGData['SGCen'])*(int(SGData['SGInv'])+1)/Jdup
747         
748    return GetKNsym(str(Isym)),Mult
749   
750def ElemPosition(SGData):
751    ''' Under development
752    Object here is to return a list of symmetry element types and locations suitable
753    for say drawing them.
754    So far I have the element type... getting all possible locations without lookup may be impossible!
755    '''
756    SymElements = []
757    Inv = SGData['SGInv']
758    Cen = SGData['SGCen']
759    eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
760    # get operators & expand if centrosymmetric
761    Ops = SGData['SGOps']
762    opM = np.array([op[0].T for op in Ops])
763    opT = np.array([op[1] for op in Ops])
764    if Inv:
765        opM = np.concatenate((opM,-opM))
766        opT = np.concatenate((opT,-opT))
767    opMT = zip(opM,opT)
768    for M,T in opMT[1:]:        #skip I
769        Dt = int(nl.det(M))
770        Tr = int(np.trace(M))
771        Dt = -(Dt-1)/2
772        Es = eleSym[Tr][Dt]
773        if Dt:              #rotation-inversion
774            I = np.eye(3)
775            if Tr == 1:     #mirrors/glides
776                if np.any(T):       #glide
777                    M2 = np.inner(M,M)
778                    MT = np.inner(M,T)+T
779                    print 'glide',Es,MT
780                    print M2
781                else:               #mirror
782                    print 'mirror',Es,T
783                    print I-M
784                X = [-1,-1,-1]
785            elif Tr == -3:  # pure inversion
786                X = np.inner(nl.inv(I-M),T)
787                print 'inversion',Es,X
788            else:           #other rotation-inversion
789                M2 = np.inner(M,M)
790                MT = np.inner(M,T)+T
791                print 'rot-inv',Es,MT
792                print M2
793                X = [-1,-1,-1]
794               
795           
796           
797        else:               #rotations
798            print 'rotation',Es
799            X = [-1,-1,-1]
800        #SymElements.append([Es,X])
801       
802    return #SymElements
803   
804def ApplyStringOps(A,SGData,X,Uij=[]):
805    SGOps = SGData['SGOps']
806    SGCen = SGData['SGCen']
807    Ax = A.split('+')
808    Ax[0] = int(Ax[0])
809    iC = 0
810    if Ax[0] < 0:
811        iC = 1
812    Ax[0] = abs(Ax[0])
813    nA = Ax[0]%100-1
814    cA = Ax[0]/100
815    Cen = SGCen[cA]
816    M,T = SGOps[nA]
817    if len(Ax)>1:
818        cellA = Ax[1].split(',')
819        cellA = np.array([int(a) for a in cellA])
820    else:
821        cellA = np.zeros(3)
822    newX = (1-2*iC)*(Cen+np.inner(M,X)+T)+cellA
823    if len(Uij):
824        U = Uij2U(Uij)
825        U = np.inner(M,np.inner(U,M).T)
826        newUij = U2Uij(U)
827        return [newX,newUij]
828    else:
829        return newX
830       
831def StringOpsProd(A,B,SGData):
832    ''' Find A*B where A & B are in strings '-' + '100*c+n' + '+ijk'
833    where '-' indicates inversion, c(>0) is the cell centering operator,
834    n is operator number from SgOps and ijk are unit cell translations (each may be <0).
835    Should return resultant string - C.
836        SGData - dictionary using entries:
837             'SGCen': cell centering vectors [0,0,0] at least
838             'SGOps': symmetry operations as [M,T] so that M*x+T = x'
839    '''
840    SGOps = SGData['SGOps']
841    SGCen = SGData['SGCen']
842    #1st split out the cell translation part & work on the operator parts
843    Ax = A.split('+'); Bx = B.split('+')
844    Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
845    iC = 0
846    if Ax[0]*Bx[0] < 0:
847        iC = 1
848    Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
849    nA = Ax[0]%100-1;  nB = Bx[0]%100-1
850    cA = Ax[0]/100;  cB = Bx[0]/100
851    Cen = (SGCen[cA]+SGCen[cB])%1.0
852    cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
853    Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
854    Mc = np.inner(Ma,Mb.T)
855#    print Ma,Mb,Mc
856    Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
857#    print Ta,Tb,Tc
858#    print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
859    nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
860    #now the cell translation part
861    if len(Ax)>1:
862        cellA = Ax[1].split(',')
863        cellA = [int(a) for a in cellA]
864    else:
865        cellA = [0,0,0]
866    if len(Bx)>1:
867        cellB = Bx[1].split(',')
868        cellB = [int(b) for b in cellB]
869    else:
870        cellB = [0,0,0]
871    cellC = np.add(cellA,cellB)
872    C = str(((nC+1)+(100*cC))*(1-2*iC))+'+'+ \
873        str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
874    return C
875           
876def U2Uij(U):
877    #returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
878    return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
879   
880def Uij2U(Uij):
881    #returns the thermal motion tensor U from Uij as numpy array
882    return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
883   
884'''A list of space groups as ordered and named in the pre-2002 International
885Tables Volume A, except that spaces are used following the GSAS convention to
886separate the different crystallographic directions.
887Note that the symmetry codes here will recognize many non-standard space group
888symbols with different settings.
889'''
890spglist = {
891    'triclinic' : ('P 1','P -1',), # 1-2
892    'monoclinic': ('P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
893                   'C 2/m','P 2/c','P 21/c','C 2/c',), #3-15
894    'orthorhombic': ('P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21','C 2 2 21',
895                     'C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
896                     'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
897                     'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
898                     'C m m 2','C m c 21','C c c 2','A m m 2','A b m 2',
899                     'A m a 2','A b a 2','F m m 2','F d d 2','I m m 2',
900                     'I b a 2','I m a 2','P m m m','P n n n','P c c m',
901                     'P b a n','P m m a','P n n a','P m n a','P c c a',
902                     'P b a m','P c c n','P b c m','P n n m','P m m n',
903                     'P b c n','P b c a','P n m a','C m c m','C m c a',
904                     'C m m m','C c c m','C m m a','C c c a','F m m m',
905                     'F d d d','I m m m','I b a m','I b c a','I m m a',), #16-74
906    'tetragonal': ('P 4','P 41','P 42','P 43','I 4','I 41','P -4','I -4',
907                   'P 4/m','P 42/m','P 4/n','P 42/n','I 4/m','I 41/a',
908                   'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
909                   'P 42 21 2','P 43 2 2','P 43 21 2','I 4 2 2','I 41 2 2',
910                   'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c',
911                   'P 4 n c','P 42 m c','P 42 b c','I 4 m m','I 4 c m',
912                   'I 41 m d','I 41 c d','P -4 2 m','P -4 2 c','P -4 21 m',
913                   'P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2','P -4 n 2',
914                   'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d','P 4/m m m',
915                   'P 4/m c c','P 4/n b m','P 4/n n c','P 4/m b m','P 4/m n c',
916                   'P 4/n m m','P 4/n c c','P 42/m m c','P 42/m c m',
917                   'P 42/n b c','P 42/n n m','P 42/m b c','P 42/m n m',
918                   'P 42/n m c','P 42/n c m','I 4/m m m','I 4/m c m',
919                   'I 41/a m d','I 41/a c d',), # 75-142
920    'trigonal': ('P 3','P 31','P 32','R 3','P -3','R -3','P 3 1 2','P 3 2 1',
921                 'P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1','R 3 2', 'P 3 m 1',
922                 'P 3 1 m','P 3 c 1','P 3 1 c','R 3 m','R 3 c','P -3 1 m',
923                 'P -3 1 c','P -3 m 1','P -3 c 1','R -3 m','R -3 c',), #143-167
924    'hexagonal': ('P 6','P 61','P 65','P 62','P 64','P 63','P -6','P 6/m',
925                  'P 63/m','P 6 2 2','P 61 2 2','P 65 2 2','P 62 2 2',
926                  'P 64 2 2','P 63 2 2','P 6 m m','P 6 c c','P 63 c m',
927                  'P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m','P -6 2 c',
928                  'P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',), #144-194
929    'cubic': ('P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m -3','P n -3',
930              'F m -3','F d -3','I m -3','P a -3','I a -3','P 4 3 2','P 42 3 2',
931              'F 4 3 2','F 41 3 2','I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2',
932              'P -4 3 m','F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
933              'P m -3 m','P n -3 n','P m -3 n','P n -3 m','F m -3 m','F m -3 c',
934              'F d -3 m','F d -3 c','I m -3 m','I a -3 d',), #195-230
935}
936'A few non-standard space groups for test use'
937nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h', 
938                      'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
939                      'R 3 c r','R -3 c r','R -3 m r',),
940'''A list of orthorhombic space groups that were renamed in the 2002 Volume A,
941along with the pre-2002 name. The e designates a double glide-plane'''
942sgequiv_2002_orthorhombic= (('A e m 2', 'A b m 2',),
943                            ('A e a 2', 'A b a 2',),
944                            ('C m c e', 'C m c a',),
945                            ('C m m e', 'C m m a',),
946                            ('C c c e', 'C c c a'),)
947'''Use the space groups types in this order to list the symbols in the
948order they are listed in the International Tables, vol. A'''
949symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 
950               'trigonal', 'hexagonal', 'cubic')
951
952# self-test materials follow. Requires files in directory testinp
953def test0():
954    '''test #0: exercise MoveToUnitCell'''
955    msg = "MoveToUnitCell failed"
956    assert (MoveToUnitCell([1,2,3]) == [0,0,0]).all, msg
957    assert (MoveToUnitCell([2,-1,-2]) == [0,0,0]).all, msg
958    assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9) < 1e-6, msg
959    assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1) < 1e-6, msg
960
961def test1():
962    ''' test #1: SpcGroup and SGPrint against previous results'''
963    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
964    if ospath.exists(testdir):
965        if testdir not in sys.path: sys.path.insert(0,testdir)
966    import spctestinp
967    def CompareSpcGroup(spc, referr, refdict, reflist): 
968        'Compare output from GSASIIspc.SpcGroup with results from a previous run'
969        # if an error is reported, the dictionary can be ignored
970        msg0 = "CompareSpcGroup failed on space group %s" % spc
971        result = SpcGroup(spc)
972        if result[0] == referr and referr > 0: return True
973        keys = result[1].keys()
974        #print result[1]['SpGrp']
975        msg = msg0 + " in list lengths"
976        assert len(keys) == len(refdict.keys()), msg
977        for key in keys:
978            if key == 'SGOps' or  key == 'SGCen':
979                msg = msg0 + (" in key %s length" % key)
980                assert len(refdict[key]) == len(result[1][key]), msg
981                for i in range(len(refdict[key])):
982                    msg = msg0 + (" in key %s level 0" % key)
983                    assert np.allclose(result[1][key][i][0],refdict[key][i][0]), msg
984                    msg = msg0 + (" in key %s level 1" % key)
985                    assert np.allclose(result[1][key][i][1],refdict[key][i][1]), msg
986            else:
987                msg = msg0 + (" in key %s" % key)
988                assert result[1][key] == refdict[key], msg
989        msg = msg0 + (" in key %s reflist" % key)
990        #for (l1,l2) in zip(reflist, SGPrint(result[1])):
991        #    assert l2.replace('\t','').replace(' ','') == l1.replace(' ',''), 'SGPrint ' +msg
992        assert reflist == SGPrint(result[1]), 'SGPrint ' +msg
993    for spc in spctestinp.SGdat:
994        CompareSpcGroup(spc, 0, spctestinp.SGdat[spc], spctestinp.SGlist[spc] )
995
996def test2():
997    ''' test #2: SpcGroup against cctbx (sgtbx) computations'''
998    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
999    if ospath.exists(testdir):
1000        if testdir not in sys.path: sys.path.insert(0,testdir)
1001    import sgtbxtestinp
1002    def CompareWcctbx(spcname, cctbx_in, debug=0):
1003        'Compare output from GSASIIspc.SpcGroup with results from cctbx.sgtbx'
1004        cctbx = cctbx_in[:] # make copy so we don't delete from the original
1005        spc = (SpcGroup(spcname))[1]
1006        if debug: print spc['SpGrp']
1007        if debug: print spc['SGCen']
1008        latticetype = spcname.strip().upper()[0]
1009        # lattice type of R implies Hexagonal centering", fix the rhombohedral settings
1010        if latticetype == "R" and len(spc['SGCen']) == 1: latticetype = 'P'
1011        assert latticetype == spc['SGLatt'], "Failed: %s does not match Lattice: %s" % (spcname, spc['SGLatt'])
1012        onebar = [1]
1013        if spc['SGInv']: onebar.append(-1)
1014        for (op,off) in spc['SGOps']:
1015            for inv in onebar:
1016                for cen in spc['SGCen']:
1017                    noff = off + cen
1018                    noff = MoveToUnitCell(noff)
1019                    mult = tuple((op*inv).ravel().tolist())
1020                    if debug: print "\n%s: %s + %s" % (spcname,mult,noff)
1021                    for refop in cctbx:
1022                        if debug: print refop
1023                        # check the transform
1024                        if refop[:9] != mult: continue
1025                        if debug: print "mult match"
1026                        # check the translation
1027                        reftrans = list(refop[-3:])
1028                        reftrans = MoveToUnitCell(reftrans)
1029                        if all(abs(noff - reftrans) < 1.e-5):
1030                            cctbx.remove(refop)
1031                            break
1032                    else:
1033                        assert False, "failed on %s:\n\t %s + %s" % (spcname,mult,noff)
1034    for key in sgtbxtestinp.sgtbx:
1035        CompareWcctbx(key, sgtbxtestinp.sgtbx[key])
1036
1037def test3(): 
1038    ''' test #3: exercise SytSym (includes GetOprPtrName, GenAtom, GetKNsym)
1039     for selected space groups against info in IT Volume A '''
1040    def ExerciseSiteSym (spc, crdlist):
1041        'compare site symmetries and multiplicities for a specified space group'
1042        msg = "failed on site sym test for %s" % spc
1043        (E,S) = SpcGroup(spc)
1044        assert not E, msg
1045        for t in crdlist:
1046            symb, m = SytSym(t[0],S)
1047            if symb.strip() != t[2].strip() or m != t[1]:
1048                print spc,t[0],m,symb
1049            assert m == t[1]
1050            #assert symb.strip() == t[2].strip()
1051
1052    ExerciseSiteSym('p 1',[
1053            ((0.13,0.22,0.31),1,'1'),
1054            ((0,0,0),1,'1'),
1055            ])
1056    ExerciseSiteSym('p -1',[
1057            ((0.13,0.22,0.31),2,'1'),
1058            ((0,0.5,0),1,'-1'),
1059            ])
1060    ExerciseSiteSym('C 2/c',[
1061            ((0.13,0.22,0.31),8,'1'),
1062            ((0.0,.31,0.25),4,'2(010)'),
1063            ((0.25,.25,0.5),4,'-1'),
1064            ((0,0.5,0),4,'-1'),
1065            ])
1066    ExerciseSiteSym('p 2 2 2',[
1067            ((0.13,0.22,0.31),4,'1'),
1068            ((0,0.5,.31),2,'2(001)'),
1069            ((0.5,.31,0.5),2,'2(010)'),
1070            ((.11,0,0),2,'2(100)'),
1071            ((0,0.5,0),1,'222'),
1072            ])
1073    ExerciseSiteSym('p 4/n',[
1074            ((0.13,0.22,0.31),8,'1'),
1075            ((0.25,0.75,.31),4,'2(001)'),
1076            ((0.5,0.5,0.5),4,'-1'),
1077            ((0,0.5,0),4,'-1'),
1078            ((0.25,0.25,.31),2,'4(001)'),
1079            ((0.25,.75,0.5),2,'-4(001)'),
1080            ((0.25,.75,0.0),2,'-4(001)'),
1081            ])
1082    ExerciseSiteSym('p 31 2 1',[
1083            ((0.13,0.22,0.31),6,'1'),
1084            ((0.13,0.0,0.833333333),3,'2(100)'),
1085            ((0.13,0.13,0.),3,'2(110)'),
1086            ])
1087    ExerciseSiteSym('R 3 c',[
1088            ((0.13,0.22,0.31),18,'1'),
1089            ((0.0,0.0,0.31),6,'3'),
1090            ])
1091    ExerciseSiteSym('R 3 c R',[
1092            ((0.13,0.22,0.31),6,'1'),
1093            ((0.31,0.31,0.31),2,'3(111)'),
1094            ])
1095    ExerciseSiteSym('P 63 m c',[
1096            ((0.13,0.22,0.31),12,'1'),
1097            ((0.11,0.22,0.31),6,'m(100)'),
1098            ((0.333333,0.6666667,0.31),2,'3m(100)'),
1099            ((0,0,0.31),2,'3m(100)'),
1100            ])
1101    ExerciseSiteSym('I a -3',[
1102            ((0.13,0.22,0.31),48,'1'),
1103            ((0.11,0,0.25),24,'2(100)'),
1104            ((0.11,0.11,0.11),16,'3(111)'),
1105            ((0,0,0),8,'-3(111)'),
1106            ])
1107
1108if __name__ == '__main__':
1109    test0()
1110    test1()
1111    test2()
1112    test3()
1113    print "OK"
Note: See TracBrowser for help on using the repository browser.