source: trunk/GSASIIspc.py @ 404

Last change on this file since 404 was 404, checked in by vondreele, 10 years ago

put error bars on seq refinement plots
add cubic elastic strain coeff
make sure results are saved from seq refinements

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