source: trunk/GSASIIstrMath.py @ 1025

Last change on this file since 1025 was 1025, checked in by vondreele, 9 years ago

fix excluded region problems

File size: 99.5 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2013-05-17 12:13:15 -0500 (Fri, 17 May 2013) $
8# $Author: vondreele $
9# $Revision: 920 $
10# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 920 2013-05-17 17:13:15Z vondreele $
12########### SVN repository information ###################
13import time
14import math
15import copy
16import numpy as np
17import numpy.ma as ma
18import numpy.linalg as nl
19import scipy.optimize as so
20import GSASIIpath
21GSASIIpath.SetVersionNumber("$Revision: 920 $")
22import GSASIIElem as G2el
23import GSASIIlattice as G2lat
24import GSASIIspc as G2spc
25import GSASIIpwd as G2pwd
26import GSASIImapvars as G2mv
27import GSASIImath as G2mth
28
29sind = lambda x: np.sin(x*np.pi/180.)
30cosd = lambda x: np.cos(x*np.pi/180.)
31tand = lambda x: np.tan(x*np.pi/180.)
32asind = lambda x: 180.*np.arcsin(x)/np.pi
33acosd = lambda x: 180.*np.arccos(x)/np.pi
34atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
35   
36ateln2 = 8.0*math.log(2.0)
37
38################################################################################
39##### Rigid Body Models
40################################################################################
41       
42def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
43    ''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
44    current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
45    '''
46    atxIds = ['Ax:','Ay:','Az:']
47    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
48    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
49    if not RBIds['Vector'] and not RBIds['Residue']:
50        return
51    VRBIds = RBIds['Vector']
52    RRBIds = RBIds['Residue']
53    if Update:
54        RBData = rigidbodyDict
55    else:
56        RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
57    if RBIds['Vector']:                       # first update the vector magnitudes
58        VRBData = RBData['Vector']
59        for i,rbId in enumerate(VRBIds):
60            if VRBData[rbId]['useCount']:
61                for j in range(len(VRBData[rbId]['VectMag'])):
62                    name = '::RBV;'+str(j)+':'+str(i)
63                    VRBData[rbId]['VectMag'][j] = parmDict[name]
64    for phase in Phases:
65        Phase = Phases[phase]
66        General = Phase['General']
67        cell = General['Cell'][1:7]
68        Amat,Bmat = G2lat.cell2AB(cell)
69        AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
70        pfx = str(Phase['pId'])+'::'
71        if Update:
72            RBModels = Phase['RBModels']
73        else:
74            RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
75        for irb,RBObj in enumerate(RBModels.get('Vector',[])):
76            jrb = VRBIds.index(RBObj['RBId'])
77            rbsx = str(irb)+':'+str(jrb)
78            for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
79                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
80            for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
81                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
82            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
83            TLS = RBObj['ThermalMotion']
84            if 'T' in TLS[0]:
85                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
86                    TLS[1][i] = parmDict[pfx+pt+rbsx]
87            if 'L' in TLS[0]:
88                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
89                    TLS[1][i+6] = parmDict[pfx+pt+rbsx]
90            if 'S' in TLS[0]:
91                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
92                    TLS[1][i+12] = parmDict[pfx+pt+rbsx]
93            if 'U' in TLS[0]:
94                TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
95            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
96            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
97            for i,x in enumerate(XYZ):
98                atId = RBObj['Ids'][i]
99                for j in [0,1,2]:
100                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
101                if UIJ[i][0] == 'A':
102                    for j in range(6):
103                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
104                elif UIJ[i][0] == 'I':
105                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
106           
107        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
108            jrb = RRBIds.index(RBObj['RBId'])
109            rbsx = str(irb)+':'+str(jrb)
110            for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
111                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
112            for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
113                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]               
114            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
115            TLS = RBObj['ThermalMotion']
116            if 'T' in TLS[0]:
117                for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
118                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
119            if 'L' in TLS[0]:
120                for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
121                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
122            if 'S' in TLS[0]:
123                for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
124                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
125            if 'U' in TLS[0]:
126                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
127            for itors,tors in enumerate(RBObj['Torsions']):
128                tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
129            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
130            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
131            for i,x in enumerate(XYZ):
132                atId = RBObj['Ids'][i]
133                for j in [0,1,2]:
134                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
135                if UIJ[i][0] == 'A':
136                    for j in range(6):
137                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
138                elif UIJ[i][0] == 'I':
139                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
140                   
141def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
142    'Needs a doc string'
143    atxIds = ['dAx:','dAy:','dAz:']
144    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
145    TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
146    LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
147    SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
148    PIds = ['Px:','Py:','Pz:']
149    OIds = ['Oa:','Oi:','Oj:','Ok:']
150    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
151    if not RBIds['Vector'] and not RBIds['Residue']:
152        return
153    VRBIds = RBIds['Vector']
154    RRBIds = RBIds['Residue']
155    RBData = rigidbodyDict
156    for item in parmDict:
157        if 'RB' in item:
158            dFdvDict[item] = 0.        #NB: this is a vector which is no. refl. long & must be filled!
159    General = Phase['General']
160    cell = General['Cell'][1:7]
161    Amat,Bmat = G2lat.cell2AB(cell)
162    rpd = np.pi/180.
163    rpd2 = rpd**2
164    g = nl.inv(np.inner(Bmat,Bmat))
165    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
166        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
167    AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
168    pfx = str(Phase['pId'])+'::'
169    RBModels =  Phase['RBModels']
170    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
171        VModel = RBData['Vector'][RBObj['RBId']]
172        Q = RBObj['Orient'][0]
173        Pos = RBObj['Orig'][0]
174        jrb = VRBIds.index(RBObj['RBId'])
175        rbsx = str(irb)+':'+str(jrb)
176        dXdv = []
177        for iv in range(len(VModel['VectMag'])):
178            dCdv = []
179            for vec in VModel['rbVect'][iv]:
180                dCdv.append(G2mth.prodQVQ(Q,vec))
181            dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
182        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
183        for ia,atId in enumerate(RBObj['Ids']):
184            atNum = AtLookup[atId]
185            dx = 0.00001
186            for iv in range(len(VModel['VectMag'])):
187                for ix in [0,1,2]:
188                    dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
189            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
190                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
191            for iv in range(4):
192                Q[iv] -= dx
193                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
194                Q[iv] += 2.*dx
195                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
196                Q[iv] -= dx
197                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
198                for ix in [0,1,2]:
199                    dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
200            X = G2mth.prodQVQ(Q,Cart[ia])
201            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
202            dFdu = G2lat.U6toUij(dFdu.T)
203            dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))           
204            dFdu = G2lat.UijtoU6(dFdu)
205            atNum = AtLookup[atId]
206            if 'T' in RBObj['ThermalMotion'][0]:
207                for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
208                    dFdvDict[pfx+name+rbsx] += dFdu[i]
209            if 'L' in RBObj['ThermalMotion'][0]:
210                dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
211                dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
212                dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
213                dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
214                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
215                dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2*(-dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
216                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
217                dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2*(-dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
218                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
219            if 'S' in RBObj['ThermalMotion'][0]:
220                dFdvDict[pfx+'RBVS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
221                dFdvDict[pfx+'RBVS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
222                dFdvDict[pfx+'RBVS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
223                dFdvDict[pfx+'RBVS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
224                dFdvDict[pfx+'RBVS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
225                dFdvDict[pfx+'RBVS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
226                dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
227                dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
228            if 'U' in RBObj['ThermalMotion'][0]:
229                dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
230
231
232    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
233        Q = RBObj['Orient'][0]
234        Pos = RBObj['Orig'][0]
235        jrb = RRBIds.index(RBObj['RBId'])
236        torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
237        rbsx = str(irb)+':'+str(jrb)
238        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
239        for itors,tors in enumerate(RBObj['Torsions']):     #derivative error?
240            tname = pfx+'RBRTr;'+str(itors)+':'+rbsx           
241            orId,pvId = torData[itors][:2]
242            pivotVec = Cart[orId]-Cart[pvId]
243            QA = G2mth.AVdeg2Q(-0.001,pivotVec)
244            QB = G2mth.AVdeg2Q(0.001,pivotVec)
245            for ir in torData[itors][3]:
246                atNum = AtLookup[RBObj['Ids'][ir]]
247                rVec = Cart[ir]-Cart[pvId]
248                dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
249                dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
250                for ix in [0,1,2]:
251                    dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
252        for ia,atId in enumerate(RBObj['Ids']):
253            atNum = AtLookup[atId]
254            dx = 0.00001
255            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
256                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
257            for iv in range(4):
258                Q[iv] -= dx
259                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
260                Q[iv] += 2.*dx
261                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
262                Q[iv] -= dx
263                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
264                for ix in [0,1,2]:
265                    dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
266            X = G2mth.prodQVQ(Q,Cart[ia])
267            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
268            dFdu = G2lat.U6toUij(dFdu.T)
269            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
270            dFdu = G2lat.UijtoU6(dFdu)
271            atNum = AtLookup[atId]
272            if 'T' in RBObj['ThermalMotion'][0]:
273                for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
274                    dFdvDict[pfx+name+rbsx] += dFdu[i]
275            if 'L' in RBObj['ThermalMotion'][0]:
276                dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
277                dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
278                dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
279                dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
280                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
281                dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2*(dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
282                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
283                dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2*(dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
284                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
285            if 'S' in RBObj['ThermalMotion'][0]:
286                dFdvDict[pfx+'RBRS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
287                dFdvDict[pfx+'RBRS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
288                dFdvDict[pfx+'RBRS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
289                dFdvDict[pfx+'RBRS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
290                dFdvDict[pfx+'RBRS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
291                dFdvDict[pfx+'RBRS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
292                dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
293                dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
294            if 'U' in RBObj['ThermalMotion'][0]:
295                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
296   
297################################################################################
298##### Penalty & restraint functions
299################################################################################
300
301def penaltyFxn(HistoPhases,parmDict,varyList):
302    'Needs a doc string'
303    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
304    pNames = []
305    pVals = []
306    pWt = []
307    negWt = {}
308    for phase in Phases:
309        pId = Phases[phase]['pId']
310        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
311        General = Phases[phase]['General']
312        textureData = General['SH Texture']
313        SGData = General['SGData']
314        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
315        cell = General['Cell'][1:7]
316        Amat,Bmat = G2lat.cell2AB(cell)
317        if phase not in restraintDict:
318            continue
319        phaseRest = restraintDict[phase]
320        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
321            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
322            ['ChemComp','Sites'],['Texture','HKLs']]
323        for name,rest in names:
324            itemRest = phaseRest[name]
325            if itemRest[rest] and itemRest['Use']:
326                wt = itemRest['wtFactor']
327                if name in ['Bond','Angle','Plane','Chiral']:
328                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
329                        pNames.append(str(pId)+':'+name+':'+str(i))
330                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
331                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
332                        if name == 'Bond':
333                            calc = G2mth.getRestDist(XYZ,Amat)
334                        elif name == 'Angle':
335                            calc = G2mth.getRestAngle(XYZ,Amat)
336                        elif name == 'Plane':
337                            calc = G2mth.getRestPlane(XYZ,Amat)
338                        elif name == 'Chiral':
339                            calc = G2mth.getRestChiral(XYZ,Amat)
340                        pVals.append(obs-calc)
341                        pWt.append(wt/esd**2)
342                elif name in ['Torsion','Rama']:
343                    coeffDict = itemRest['Coeff']
344                    for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
345                        pNames.append(str(pId)+':'+name+':'+str(i))
346                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
347                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
348                        if name == 'Torsion':
349                            tor = G2mth.getRestTorsion(XYZ,Amat)
350                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
351                        else:
352                            phi,psi = G2mth.getRestRama(XYZ,Amat)
353                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
354                        pVals.append(obs-calc)
355                        pWt.append(wt/esd**2)
356                elif name == 'ChemComp':
357                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
358                        pNames.append(str(pId)+':'+name+':'+str(i))
359                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
360                        frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
361                        calc = np.sum(mul*frac*factors)
362                        pVals.append(obs-calc)
363                        pWt.append(wt/esd**2)                   
364                elif name == 'Texture':
365                    SHkeys = textureData['SH Coeff'][1].keys()
366                    SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
367                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
368                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
369                    for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
370                        PH = np.array(hkl)
371                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
372                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
373                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
374                        Z1 = -ma.masked_greater(Z,0.0)
375                        IndZ1 = np.array(ma.nonzero(Z1))
376                        for ind in IndZ1.T:
377                            pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
378                            pVals.append(Z1[ind[0]][ind[1]])
379                            pWt.append(wt/esd1**2)
380                        if ifesd2:
381                            Z2 = 1.-Z
382                            for ind in np.ndindex(grid,grid):
383                                pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
384                                pVals.append(Z1[ind[0]][ind[1]])
385                                pWt.append(wt/esd2**2)
386         
387    for item in varyList:
388        if 'PWLref' in item and parmDict[item] < 0.:
389            pId = int(item.split(':')[0])
390            if negWt[pId]:
391                pNames.append(item)
392                pVals.append(-parmDict[item])
393                pWt.append(negWt[pId])
394    pVals = np.array(pVals)
395    pWt = np.array(pWt)         #should this be np.sqrt?
396    return pNames,pVals,pWt
397   
398def penaltyDeriv(pNames,pVal,HistoPhases,parmDict,varyList):
399    'Needs a doc string'
400    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
401    pDerv = np.zeros((len(varyList),len(pVal)))
402    for phase in Phases:
403#        if phase not in restraintDict:
404#            continue
405        pId = Phases[phase]['pId']
406        General = Phases[phase]['General']
407        SGData = General['SGData']
408        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
409        cell = General['Cell'][1:7]
410        Amat,Bmat = G2lat.cell2AB(cell)
411        textureData = General['SH Texture']
412
413        SHkeys = textureData['SH Coeff'][1].keys()
414        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
415        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
416        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
417        sam = SamSym[textureData['Model']]
418        phaseRest = restraintDict.get(phase,{})
419        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
420            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
421            'ChemComp':'Sites','Texture':'HKLs'}
422        lasthkl = np.array([0,0,0])
423        for ip,pName in enumerate(pNames):
424            pnames = pName.split(':')
425            if pId == int(pnames[0]):
426                name = pnames[1]
427                if 'PWL' in pName:
428                    pDerv[varyList.index(pName)][ip] += 1.
429                    continue
430                id = int(pnames[2]) 
431                itemRest = phaseRest[name]
432                if name in ['Bond','Angle','Plane','Chiral']:
433                    indx,ops,obs,esd = itemRest[names[name]][id]
434                    dNames = []
435                    for ind in indx:
436                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
437                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
438                    if name == 'Bond':
439                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
440                    elif name == 'Angle':
441                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
442                    elif name == 'Plane':
443                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
444                    elif name == 'Chiral':
445                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
446                elif name in ['Torsion','Rama']:
447                    coffDict = itemRest['Coeff']
448                    indx,ops,cofName,esd = itemRest[names[name]][id]
449                    dNames = []
450                    for ind in indx:
451                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
452                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
453                    if name == 'Torsion':
454                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
455                    else:
456                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
457                elif name == 'ChemComp':
458                    indx,factors,obs,esd = itemRest[names[name]][id]
459                    dNames = []
460                    for ind in indx:
461                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
462                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
463                        deriv = mul*factors
464                elif 'Texture' in name:
465                    deriv = []
466                    dNames = []
467                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
468                    hkl = np.array(hkl)
469                    if np.any(lasthkl-hkl):
470                        PH = np.array(hkl)
471                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
472                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
473                        lasthkl = copy.copy(hkl)                       
474                    if 'unit' in name:
475                        pass
476                    else:
477                        gam = float(pnames[3])
478                        psi = float(pnames[4])
479                        for SHname in ODFln:
480                            l,m,n = eval(SHname[1:])
481                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
482                            dNames += [str(pId)+'::'+SHname]
483                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
484                for dName,drv in zip(dNames,deriv):
485                    try:
486                        ind = varyList.index(dName)
487                        pDerv[ind][ip] += drv
488                    except ValueError:
489                        pass
490    return pDerv
491
492################################################################################
493##### Function & derivative calculations
494################################################################################       
495                   
496def GetAtomFXU(pfx,calcControls,parmDict):
497    'Needs a doc string'
498    Natoms = calcControls['Natoms'][pfx]
499    Tdata = Natoms*[' ',]
500    Mdata = np.zeros(Natoms)
501    IAdata = Natoms*[' ',]
502    Fdata = np.zeros(Natoms)
503    FFdata = []
504    BLdata = []
505    Xdata = np.zeros((3,Natoms))
506    dXdata = np.zeros((3,Natoms))
507    Uisodata = np.zeros(Natoms)
508    Uijdata = np.zeros((6,Natoms))
509    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
510        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
511        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
512        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
513        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
514    for iatm in range(Natoms):
515        for key in keys:
516            parm = pfx+key+str(iatm)
517            if parm in parmDict:
518                keys[key][iatm] = parmDict[parm]
519    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
520   
521def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
522    ''' Compute structure factors for all h,k,l for phase
523    puts the result, F^2, in each ref[8] in refList
524    input:
525   
526    :param list refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
527    :param np.array G:      reciprocal metric tensor
528    :param str pfx:    phase id string
529    :param dict SGData: space group info. dictionary output from SpcGroup
530    :param dict calcControls:
531    :param dict ParmDict:
532
533    '''       
534    twopi = 2.0*np.pi
535    twopisq = 2.0*np.pi**2
536    phfx = pfx.split(':')[0]+hfx
537    ast = np.sqrt(np.diag(G))
538    Mast = twopisq*np.multiply.outer(ast,ast)
539    FFtables = calcControls['FFtables']
540    BLtables = calcControls['BLtables']
541    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
542    FF = np.zeros(len(Tdata))
543    if 'N' in parmDict[hfx+'Type']:
544        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
545    else:
546        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
547        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
548    Uij = np.array(G2lat.U6toUij(Uijdata))
549    bij = Mast*Uij.T
550    for refl in refList:
551        fbs = np.array([0,0])
552        H = refl[:3]
553        SQ = 1./(2.*refl[4])**2
554        SQfactor = 4.0*SQ*twopisq
555        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
556        if not len(refl[-1]):                #no form factors
557            if 'N' in parmDict[hfx+'Type']:
558                refl[-1] = G2el.getBLvalues(BLtables)
559            else:       #'X'
560                refl[-1] = G2el.getFFvalues(FFtables,SQ)
561        for i,El in enumerate(Tdata):
562            FF[i] = refl[-1][El]           
563        Uniq = refl[11]
564        phi = refl[12]
565        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
566        sinp = np.sin(phase)
567        cosp = np.cos(phase)
568        occ = Mdata*Fdata/len(Uniq)
569        biso = -SQfactor*Uisodata
570        Tiso = np.where(biso<1.,np.exp(biso),1.0)
571        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
572        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
573        Tcorr = Tiso*Tuij
574        fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
575        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
576        if not SGData['SGInv']:
577            fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
578            fbs = np.sum(np.sum(fb,axis=1),axis=1)
579        fasq = fas**2
580        fbsq = fbs**2        #imaginary
581        refl[9] = np.sum(fasq)+np.sum(fbsq)
582        refl[10] = atan2d(fbs[0],fas[0])
583   
584def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
585    'Needs a doc string'
586    twopi = 2.0*np.pi
587    twopisq = 2.0*np.pi**2
588    phfx = pfx.split(':')[0]+hfx
589    ast = np.sqrt(np.diag(G))
590    Mast = twopisq*np.multiply.outer(ast,ast)
591    FFtables = calcControls['FFtables']
592    BLtables = calcControls['BLtables']
593    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
594    FF = np.zeros(len(Tdata))
595    if 'N' in parmDict[hfx+'Type']:
596        FP = 0.
597        FPP = 0.
598    else:
599        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
600        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
601    Uij = np.array(G2lat.U6toUij(Uijdata))
602    bij = Mast*Uij.T
603    dFdvDict = {}
604    dFdfr = np.zeros((len(refList),len(Mdata)))
605    dFdx = np.zeros((len(refList),len(Mdata),3))
606    dFdui = np.zeros((len(refList),len(Mdata)))
607    dFdua = np.zeros((len(refList),len(Mdata),6))
608    dFdbab = np.zeros((len(refList),2))
609    for iref,refl in enumerate(refList):
610        H = np.array(refl[:3])
611        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
612        SQfactor = 8.0*SQ*np.pi**2
613        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
614        Bab = parmDict[phfx+'BabA']*dBabdA
615        for i,El in enumerate(Tdata):           
616            FF[i] = refl[-1][El]           
617        Uniq = refl[11]
618        phi = refl[12]
619        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
620        sinp = np.sin(phase)
621        cosp = np.cos(phase)
622        occ = Mdata*Fdata/len(Uniq)
623        biso = -SQfactor*Uisodata
624        Tiso = np.where(biso<1.,np.exp(biso),1.0)
625        HbH = -np.inner(H,np.inner(bij,H))
626        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
627        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
628        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
629        Tcorr = Tiso*Tuij
630        fot = (FF+FP-Bab)*occ*Tcorr
631        fotp = FPP*occ*Tcorr
632        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
633        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
634       
635        fas = np.sum(np.sum(fa,axis=1),axis=1)
636        fbs = np.sum(np.sum(fb,axis=1),axis=1)
637        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
638        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
639        #sum below is over Uniq
640        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
641        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
642        dfadui = np.sum(-SQfactor*fa,axis=2)
643        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
644        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
645        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
646        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
647        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
648        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
649        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
650        dFdbab[iref] = np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
651        if not SGData['SGInv']:
652            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
653            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
654            dfbdui = np.sum(-SQfactor*fb,axis=2)
655            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
656            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
657            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
658            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
659            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
660            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
661            dFdbab[iref] += np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
662        #loop over atoms - each dict entry is list of derivatives for all the reflections
663    for i in range(len(Mdata)):     
664        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
665        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
666        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
667        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
668        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
669        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
670        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
671        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
672        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
673        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
674        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
675        dFdvDict[pfx+'BabA'] = dFdbab.T[0]
676        dFdvDict[pfx+'BabU'] = dFdbab.T[1]
677    return dFdvDict
678   
679def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
680    ''' Single crystal extinction function; puts correction in ref[13] and returns
681    corrections needed for derivatives
682    '''
683    ref[13] = 1.0
684    dervCor = 1.0
685    dervDict = {}
686    if calcControls[phfx+'EType'] != 'None':
687        cos2T = 1.0-0.5*(parmDict[hfx+'Lam']/ref[4])**2         #cos(2theta)
688        if 'SXC' in parmDict[hfx+'Type']:
689            AV = 7.9406e5/parmDict[pfx+'Vol']**2
690            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
691            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
692        elif 'SNT' in parmDict[hfx+'Type']:
693            AV = 1.e7/parmDict[pfx+'Vol']**2
694            PL = 1./(4.*refl[4]**2)
695            P12 = 1.0
696        elif 'SNC' in parmDict[hfx+'Type']:
697            AV = 1.e7/parmDict[pfx+'Vol']**2
698            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
699            P12 = 1.0
700           
701        PLZ = AV*P12*parmDict[hfx+'Lam']**2*ref[7]
702        if 'Primary' in calcControls[phfx+'EType']:
703            PLZ *= 1.5
704        else:
705            PLZ *= calcControls[phfx+'Tbar']
706                       
707        if 'Primary' in calcControls[phfx+'EType']:
708            PSIG = parmDict[phfx+'Ep']
709        elif 'I & II' in calcControls[phfx+'EType']:
710            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
711        elif 'Type II' in calcControls[phfx+'EType']:
712            PSIG = parmDict[phfx+'Es']
713        else:       # 'Secondary Type I'
714            PSIG = parmDict[phfx+'Eg']/PL
715           
716        AG = 0.58+0.48*cos2T+0.24*cos2T**2
717        AL = 0.025+0.285*cos2T
718        BG = 0.02-0.025*cos2T
719        BL = 0.15-0.2*(0.75-cos2T)**2
720        if cos2T < 0.:
721            BL = -0.45*cos2T
722        CG = 2.
723        CL = 2.
724        PF = PLZ*PSIG
725       
726        if 'Gaussian' in calcControls[phfx+'EApprox']:
727            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
728            extCor = np.sqrt(PF4)
729            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
730        else:
731            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
732            extCor = np.sqrt(PF4)
733            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
734
735        dervCor = (1.+PF)*PF3
736        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
737            dervDict[phfx+'Ep'] = -ref[7]*PLZ*PF3
738        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
739            dervDict[phfx+'Es'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
740        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
741            dervDict[phfx+'Eg'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
742               
743        ref[13] = 1./extCor
744    return dervCor,dervDict
745       
746   
747def Dict2Values(parmdict, varylist):
748    '''Use before call to leastsq to setup list of values for the parameters
749    in parmdict, as selected by key in varylist'''
750    return [parmdict[key] for key in varylist] 
751   
752def Values2Dict(parmdict, varylist, values):
753    ''' Use after call to leastsq to update the parameter dictionary with
754    values corresponding to keys in varylist'''
755    parmdict.update(zip(varylist,values))
756   
757def GetNewCellParms(parmDict,varyList):
758    'Needs a doc string'
759    newCellDict = {}
760    Anames = ['A'+str(i) for i in range(6)]
761    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
762    for item in varyList:
763        keys = item.split(':')
764        if keys[2] in Ddict:
765            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
766            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
767            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
768    return newCellDict          # is e.g. {'0::D11':A0+D11}
769   
770def ApplyXYZshifts(parmDict,varyList):
771    '''
772    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
773   
774    :param dict parmDict: parameter dictionary
775    :param list varyList: list of variables
776    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
777
778    '''
779    newAtomDict = {}
780    for item in parmDict:
781        if 'dA' in item:
782            parm = ''.join(item.split('d'))
783            parmDict[parm] += parmDict[item]
784            newAtomDict[item] = [parm,parmDict[parm]]
785    return newAtomDict
786   
787def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
788    'Spherical harmonics texture'
789    IFCoup = 'Bragg' in calcControls[hfx+'instType']
790    odfCor = 1.0
791    H = refl[:3]
792    cell = G2lat.Gmat2cell(g)
793    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
794    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
795    phi,beta = G2lat.CrsAng(H,cell,SGData)
796    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
797    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
798    for item in SHnames:
799        L,M,N = eval(item.strip('C'))
800        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
801        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
802        Lnorm = G2lat.Lnorm(L)
803        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
804    return odfCor
805   
806def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
807    'Spherical harmonics texture derivatives'
808    FORPI = 4.0*np.pi
809    IFCoup = 'Bragg' in calcControls[hfx+'instType']
810    odfCor = 1.0
811    dFdODF = {}
812    dFdSA = [0,0,0]
813    H = refl[:3]
814    cell = G2lat.Gmat2cell(g)
815    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
816    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
817    phi,beta = G2lat.CrsAng(H,cell,SGData)
818    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
819    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
820    for item in SHnames:
821        L,M,N = eval(item.strip('C'))
822        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
823        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
824        Lnorm = G2lat.Lnorm(L)
825        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
826        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
827        for i in range(3):
828            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
829    return odfCor,dFdODF,dFdSA
830   
831def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
832    'spherical harmonics preferred orientation (cylindrical symmetry only)'
833    odfCor = 1.0
834    H = refl[:3]
835    cell = G2lat.Gmat2cell(g)
836    Sangl = [0.,0.,0.]
837    if 'Bragg' in calcControls[hfx+'instType']:
838        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
839        IFCoup = True
840    else:
841        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
842        IFCoup = False
843    phi,beta = G2lat.CrsAng(H,cell,SGData)
844    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
845    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
846    for item in SHnames:
847        L,N = eval(item.strip('C'))
848        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
849        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
850    return odfCor
851   
852def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
853    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
854    FORPI = 12.5663706143592
855    odfCor = 1.0
856    dFdODF = {}
857    H = refl[:3]
858    cell = G2lat.Gmat2cell(g)
859    Sangl = [0.,0.,0.]
860    if 'Bragg' in calcControls[hfx+'instType']:
861        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
862        IFCoup = True
863    else:
864        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
865        IFCoup = False
866    phi,beta = G2lat.CrsAng(H,cell,SGData)
867    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
868    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
869    for item in SHnames:
870        L,N = eval(item.strip('C'))
871        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
872        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
873        dFdODF[phfx+item] = Kcsl*Lnorm
874    return odfCor,dFdODF
875   
876def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
877    'Needs a doc string'
878    POcorr = 1.0
879    if calcControls[phfx+'poType'] == 'MD':
880        MD = parmDict[phfx+'MD']
881        if MD != 1.0:
882            MDAxis = calcControls[phfx+'MDAxis']
883            sumMD = 0
884            for H in refl[11]:           
885                cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
886                A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
887                sumMD += A**3
888            POcorr = sumMD/len(refl[11])
889    else:   #spherical harmonics
890        if calcControls[phfx+'SHord']:
891            POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
892    return POcorr
893   
894def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
895    'Needs a doc string'
896    POcorr = 1.0
897    POderv = {}
898    if calcControls[phfx+'poType'] == 'MD':
899        MD = parmDict[phfx+'MD']
900        MDAxis = calcControls[phfx+'MDAxis']
901        sumMD = 0
902        sumdMD = 0
903        for H in refl[11]:           
904            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
905            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
906            sumMD += A**3
907            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
908        POcorr = sumMD/len(refl[11])
909        POderv[phfx+'MD'] = sumdMD/len(refl[11])
910    else:   #spherical harmonics
911        if calcControls[phfx+'SHord']:
912            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
913    return POcorr,POderv
914   
915def GetAbsorb(refl,hfx,calcControls,parmDict):
916    'Needs a doc string'
917    if 'Debye' in calcControls[hfx+'instType']:
918        return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
919    else:
920        return 1.0
921   
922def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
923    'Needs a doc string'
924    if 'Debye' in calcControls[hfx+'instType']:
925        return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
926    else:
927        return 0.0
928   
929def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
930    'Needs a doc string'
931    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
932    if 'X' in parmDict[hfx+'Type']:
933        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
934    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
935    if pfx+'SHorder' in parmDict:
936        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
937    Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
938    refl[13] = Icorr       
939   
940def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
941    'Needs a doc string'
942    dIdsh = 1./parmDict[hfx+'Scale']
943    dIdsp = 1./parmDict[phfx+'Scale']
944    if 'X' in parmDict[hfx+'Type']:
945        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
946        dIdPola /= pola
947    else:       #'N'
948        dIdPola = 0.0
949    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
950    for iPO in dIdPO:
951        dIdPO[iPO] /= POcorr
952    dFdODF = {}
953    dFdSA = [0,0,0]
954    if pfx+'SHorder' in parmDict:
955        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
956        for iSH in dFdODF:
957            dFdODF[iSH] /= odfCor
958        for i in range(3):
959            dFdSA[i] /= odfCor
960    dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
961    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
962       
963def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
964    'Needs a doc string'
965    costh = cosd(refl[5]/2.)
966    #crystallite size
967    if calcControls[phfx+'SizeType'] == 'isotropic':
968        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
969    elif calcControls[phfx+'SizeType'] == 'uniaxial':
970        H = np.array(refl[:3])
971        P = np.array(calcControls[phfx+'SizeAxis'])
972        cosP,sinP = G2lat.CosSinAngle(H,P,G)
973        Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
974        Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
975    else:           #ellipsoidal crystallites
976        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
977        H = np.array(refl[:3])
978        lenR = G2pwd.ellipseSize(H,Sij,GB)
979        Sgam = 1.8*wave/(np.pi*costh*lenR)
980    #microstrain               
981    if calcControls[phfx+'MustrainType'] == 'isotropic':
982        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
983    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
984        H = np.array(refl[:3])
985        P = np.array(calcControls[phfx+'MustrainAxis'])
986        cosP,sinP = G2lat.CosSinAngle(H,P,G)
987        Si = parmDict[phfx+'Mustrain;i']
988        Sa = parmDict[phfx+'Mustrain;a']
989        Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
990    else:       #generalized - P.W. Stephens model
991        pwrs = calcControls[phfx+'MuPwrs']
992        sum = 0
993        for i,pwr in enumerate(pwrs):
994            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
995        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
996    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
997    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
998    sig /= ateln2
999    return sig,gam
1000       
1001def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
1002    'Needs a doc string'
1003    gamDict = {}
1004    sigDict = {}
1005    costh = cosd(refl[5]/2.)
1006    tanth = tand(refl[5]/2.)
1007    #crystallite size derivatives
1008    if calcControls[phfx+'SizeType'] == 'isotropic':
1009        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1010        gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1011        sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1012    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1013        H = np.array(refl[:3])
1014        P = np.array(calcControls[phfx+'SizeAxis'])
1015        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1016        Si = parmDict[phfx+'Size;i']
1017        Sa = parmDict[phfx+'Size;a']
1018        gami = (1.8*wave/np.pi)/(Si*Sa)
1019        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1020        Sgam = gami*sqtrm
1021        gam = Sgam/costh
1022        dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si)
1023        dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa)
1024        gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1025        gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1026        sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1027        sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1028    else:           #ellipsoidal crystallites
1029        const = 1.8*wave/(np.pi*costh)
1030        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1031        H = np.array(refl[:3])
1032        lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1033        Sgam = 1.8*wave/(np.pi*costh*lenR)
1034        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1035            gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1036            sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1037    gamDict[phfx+'Size;mx'] = Sgam
1038    sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1039           
1040    #microstrain derivatives               
1041    if calcControls[phfx+'MustrainType'] == 'isotropic':
1042        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
1043        gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1044        sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1045    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1046        H = np.array(refl[:3])
1047        P = np.array(calcControls[phfx+'MustrainAxis'])
1048        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1049        Si = parmDict[phfx+'Mustrain;i']
1050        Sa = parmDict[phfx+'Mustrain;a']
1051        gami = 0.018*Si*Sa*tanth/np.pi
1052        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1053        Mgam = gami/sqtrm
1054        dsi = -gami*Si*cosP**2/sqtrm**3
1055        dsa = -gami*Sa*sinP**2/sqtrm**3
1056        gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1057        gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1058        sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1059        sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1060    else:       #generalized - P.W. Stephens model
1061        pwrs = calcControls[phfx+'MuPwrs']
1062        const = 0.018*refl[4]**2*tanth
1063        sum = 0
1064        for i,pwr in enumerate(pwrs):
1065            term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1066            sum += parmDict[phfx+'Mustrain:'+str(i)]*term
1067            gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx']
1068            sigDict[phfx+'Mustrain:'+str(i)] = \
1069                2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1070        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
1071        for i in range(len(pwrs)):
1072            sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
1073    gamDict[phfx+'Mustrain;mx'] = Mgam
1074    sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1075    return sigDict,gamDict
1076       
1077def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1078    'Needs a doc string'
1079    h,k,l = refl[:3]
1080    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1081    d = np.sqrt(dsq)
1082
1083    refl[4] = d
1084    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1085    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1086    if 'Bragg' in calcControls[hfx+'instType']:
1087        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1088            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1089    else:               #Debye-Scherrer - simple but maybe not right
1090        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1091    return pos
1092
1093def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1094    'Needs a doc string'
1095    dpr = 180./np.pi
1096    h,k,l = refl[:3]
1097    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1098    dst = np.sqrt(dstsq)
1099    pos = refl[5]-parmDict[hfx+'Zero']
1100    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1101    dpdw = const*dst
1102    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1103    dpdA *= const*wave/(2.0*dst)
1104    dpdZ = 1.0
1105    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1106    if 'Bragg' in calcControls[hfx+'instType']:
1107        dpdSh = -4.*const*cosd(pos/2.0)
1108        dpdTr = -const*sind(pos)*100.0
1109        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1110    else:               #Debye-Scherrer - simple but maybe not right
1111        dpdXd = -const*cosd(pos)
1112        dpdYd = -const*sind(pos)
1113        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1114           
1115def GetHStrainShift(refl,SGData,phfx,parmDict):
1116    'Needs a doc string'
1117    laue = SGData['SGLaue']
1118    uniq = SGData['SGUniq']
1119    h,k,l = refl[:3]
1120    if laue in ['m3','m3m']:
1121        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1122            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1123    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1124        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1125    elif laue in ['3R','3mR']:
1126        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1127    elif laue in ['4/m','4/mmm']:
1128        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1129    elif laue in ['mmm']:
1130        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1131    elif laue in ['2/m']:
1132        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1133        if uniq == 'a':
1134            Dij += parmDict[phfx+'D23']*k*l
1135        elif uniq == 'b':
1136            Dij += parmDict[phfx+'D13']*h*l
1137        elif uniq == 'c':
1138            Dij += parmDict[phfx+'D12']*h*k
1139    else:
1140        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1141            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1142    return -Dij*refl[4]**2*tand(refl[5]/2.0)
1143           
1144def GetHStrainShiftDerv(refl,SGData,phfx):
1145    'Needs a doc string'
1146    laue = SGData['SGLaue']
1147    uniq = SGData['SGUniq']
1148    h,k,l = refl[:3]
1149    if laue in ['m3','m3m']:
1150        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1151            phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1152    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1153        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1154    elif laue in ['3R','3mR']:
1155        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1156    elif laue in ['4/m','4/mmm']:
1157        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1158    elif laue in ['mmm']:
1159        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1160    elif laue in ['2/m']:
1161        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1162        if uniq == 'a':
1163            dDijDict[phfx+'D23'] = k*l
1164        elif uniq == 'b':
1165            dDijDict[phfx+'D13'] = h*l
1166        elif uniq == 'c':
1167            dDijDict[phfx+'D12'] = h*k
1168            names.append()
1169    else:
1170        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1171            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1172    for item in dDijDict:
1173        dDijDict[item] *= -refl[4]**2*tand(refl[5]/2.0)
1174    return dDijDict
1175   
1176def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1177    'Needs a doc string'
1178    histoList = Histograms.keys()
1179    histoList.sort()
1180    for histogram in histoList:
1181        if 'PWDR' in histogram[:4]:
1182            Histogram = Histograms[histogram]
1183            hId = Histogram['hId']
1184            hfx = ':%d:'%(hId)
1185            Limits = calcControls[hfx+'Limits']
1186            shl = max(parmDict[hfx+'SH/L'],0.0005)
1187            Ka2 = False
1188            kRatio = 0.0
1189            if hfx+'Lam1' in parmDict.keys():
1190                Ka2 = True
1191                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1192                kRatio = parmDict[hfx+'I(L2)/I(L1)']
1193            x,y,w,yc,yb,yd = Histogram['Data']
1194            xB = np.searchsorted(x,Limits[0])
1195            xF = np.searchsorted(x,Limits[1])
1196            ymb = np.array(y-yb)
1197            ymb = np.where(ymb,ymb,1.0)
1198            ycmb = np.array(yc-yb)
1199            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1200            refLists = Histogram['Reflection Lists']
1201            for phase in refLists:
1202                Phase = Phases[phase]
1203                pId = Phase['pId']
1204                phfx = '%d:%d:'%(pId,hId)
1205                refList = refLists[phase]
1206                sumFo = 0.0
1207                sumdF = 0.0
1208                sumFosq = 0.0
1209                sumdFsq = 0.0
1210                for refl in refList:
1211                    if 'C' in calcControls[hfx+'histType']:
1212                        yp = np.zeros_like(yb)
1213                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1214                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1215                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1216                        iFin2 = iFin
1217                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1218                        if Ka2:
1219                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1220                            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1221                            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1222                            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1223                            if not iBeg2+iFin2:       #peak below low limit - skip peak
1224                                continue
1225                            elif not iBeg2-iFin2:     #peak above high limit - done
1226                                break
1227                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
1228                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
1229                    elif 'T' in calcControls[hfx+'histType']:
1230                        print 'TOF Undefined at present'
1231                        raise Exception    #no TOF yet
1232                    Fo = np.sqrt(np.abs(refl[8]))
1233                    Fc = np.sqrt(np.abs(refl[9]))
1234                    sumFo += Fo
1235                    sumFosq += refl[8]**2
1236                    sumdF += np.abs(Fo-Fc)
1237                    sumdFsq += (refl[8]-refl[9])**2
1238                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1239                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1240                Histogram[phfx+'Nref'] = len(refList)
1241               
1242def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1243    'Needs a doc string'
1244   
1245    def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
1246        U = parmDict[hfx+'U']
1247        V = parmDict[hfx+'V']
1248        W = parmDict[hfx+'W']
1249        X = parmDict[hfx+'X']
1250        Y = parmDict[hfx+'Y']
1251        tanPos = tand(refl[5]/2.0)
1252        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
1253        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1254        sig = max(0.001,sig)
1255        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
1256        gam = max(0.001,gam)
1257        return sig,gam
1258               
1259    hId = Histogram['hId']
1260    hfx = ':%d:'%(hId)
1261    bakType = calcControls[hfx+'bakType']
1262    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1263    yc = np.zeros_like(yb)
1264       
1265    if 'C' in calcControls[hfx+'histType']:   
1266        shl = max(parmDict[hfx+'SH/L'],0.002)
1267        Ka2 = False
1268        if hfx+'Lam1' in parmDict.keys():
1269            wave = parmDict[hfx+'Lam1']
1270            Ka2 = True
1271            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1272            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1273        else:
1274            wave = parmDict[hfx+'Lam']
1275    else:
1276        print 'TOF Undefined at present'
1277        raise ValueError
1278    for phase in Histogram['Reflection Lists']:
1279        refList = Histogram['Reflection Lists'][phase]
1280        Phase = Phases[phase]
1281        pId = Phase['pId']
1282        pfx = '%d::'%(pId)
1283        phfx = '%d:%d:'%(pId,hId)
1284        hfx = ':%d:'%(hId)
1285        SGData = Phase['General']['SGData']
1286        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1287        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1288        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1289        Vst = np.sqrt(nl.det(G))    #V*
1290        if not Phase['General'].get('doPawley'):
1291            time0 = time.time()
1292            StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1293#            print 'sf calc time: %.3fs'%(time.time()-time0)
1294        time0 = time.time()
1295        for refl in refList:
1296            if 'C' in calcControls[hfx+'histType']:
1297                h,k,l = refl[:3]
1298                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1299                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
1300                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
1301                refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
1302                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
1303                refl[13] *= Vst*Lorenz
1304                if Phase['General'].get('doPawley'):
1305                    try:
1306                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1307                        refl[9] = parmDict[pInd]
1308                    except KeyError:
1309#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1310                        continue
1311                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1312                iBeg = np.searchsorted(x,refl[5]-fmin)
1313                iFin = np.searchsorted(x,refl[5]+fmax)
1314                if not iBeg+iFin:       #peak below low limit - skip peak
1315                    continue
1316                elif not iBeg-iFin:     #peak above high limit - done
1317                    break
1318                   
1319                print refl[5],ma.count(x[iBeg:iFin]),iFin-iBeg,refl[13],refl[9]
1320                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
1321                if Ka2:
1322                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1323                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1324                    iBeg = np.searchsorted(x,pos2-fmin)
1325                    iFin = np.searchsorted(x,pos2+fmax)
1326                    if not iBeg+iFin:       #peak below low limit - skip peak
1327                        continue
1328                    elif not iBeg-iFin:     #peak above high limit - done
1329                        return yc,yb
1330                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))        #and here
1331            elif 'T' in calcControls[hfx+'histType']:
1332                print 'TOF Undefined at present'
1333                raise Exception    #no TOF yet
1334#        print 'profile calc time: %.3fs'%(time.time()-time0)
1335    return yc,yb
1336   
1337def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1338    'Needs a doc string'
1339   
1340    def cellVaryDerv(pfx,SGData,dpdA): 
1341        if SGData['SGLaue'] in ['-1',]:
1342            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1343                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1344        elif SGData['SGLaue'] in ['2/m',]:
1345            if SGData['SGUniq'] == 'a':
1346                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1347            elif SGData['SGUniq'] == 'b':
1348                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1349            else:
1350                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1351        elif SGData['SGLaue'] in ['mmm',]:
1352            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1353        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1354            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1355        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1356            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1357        elif SGData['SGLaue'] in ['3R', '3mR']:
1358            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1359        elif SGData['SGLaue'] in ['m3m','m3']:
1360            return [[pfx+'A0',dpdA[0]]]
1361           
1362    # create a list of dependent variables and set up a dictionary to hold their derivatives
1363    dependentVars = G2mv.GetDependentVars()
1364    depDerivDict = {}
1365    for j in dependentVars:
1366        depDerivDict[j] = np.zeros(shape=(len(x)))
1367    #print 'dependent vars',dependentVars
1368    lenX = len(x)               
1369    hId = Histogram['hId']
1370    hfx = ':%d:'%(hId)
1371    bakType = calcControls[hfx+'bakType']
1372    dMdv = np.zeros(shape=(len(varylist),len(x)))
1373    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1374    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
1375        bBpos =varylist.index(hfx+'Back:0')
1376        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1377    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1378    for name in varylist:
1379        if 'Debye' in name:
1380            id = int(name.split(':')[-1])
1381            parm = name[:int(name.rindex(':'))]
1382            ip = names.index(parm)
1383            dMdv[varylist.index(name)] = dMddb[3*id+ip]
1384    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1385    for name in varylist:
1386        if 'BkPk' in name:
1387            id = int(name.split(':')[-1])
1388            parm = name[:int(name.rindex(':'))]
1389            ip = names.index(parm)
1390            dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1391    cw = np.diff(x)
1392    cw = np.append(cw,cw[-1])
1393    if 'C' in calcControls[hfx+'histType']:   
1394        shl = max(parmDict[hfx+'SH/L'],0.002)
1395        Ka2 = False
1396        if hfx+'Lam1' in parmDict.keys():
1397            wave = parmDict[hfx+'Lam1']
1398            Ka2 = True
1399            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1400            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1401        else:
1402            wave = parmDict[hfx+'Lam']
1403    else:
1404        print 'TOF Undefined at present'
1405        raise ValueError
1406    for phase in Histogram['Reflection Lists']:
1407        refList = Histogram['Reflection Lists'][phase]
1408        Phase = Phases[phase]
1409        SGData = Phase['General']['SGData']
1410        pId = Phase['pId']
1411        pfx = '%d::'%(pId)
1412        phfx = '%d:%d:'%(pId,hId)
1413        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1414        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1415        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1416        if not Phase['General'].get('doPawley'):
1417            time0 = time.time()
1418            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1419#            print 'sf-derv time %.3fs'%(time.time()-time0)
1420            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1421        time0 = time.time()
1422        for iref,refl in enumerate(refList):
1423            if 'C' in calcControls[hfx+'histType']:        #CW powder
1424                h,k,l = refl[:3]
1425                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1426                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1427                iBeg = np.searchsorted(x,refl[5]-fmin)
1428                iFin = np.searchsorted(x,refl[5]+fmax)
1429                if not iBeg+iFin:       #peak below low limit - skip peak
1430                    continue
1431                elif not iBeg-iFin:     #peak above high limit - done
1432                    break
1433                pos = refl[5]
1434                tanth = tand(pos/2.0)
1435                costh = cosd(pos/2.0)
1436                lenBF = iFin-iBeg
1437                dMdpk = np.zeros(shape=(6,lenBF))
1438                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
1439                for i in range(5):
1440                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[13]*refl[9]*dMdipk[i]
1441                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1442                if Ka2:
1443                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
1444                    iBeg2 = np.searchsorted(x,pos2-fmin)
1445                    iFin2 = np.searchsorted(x,pos2+fmax)
1446                    if iBeg2-iFin2:
1447                        lenBF2 = iFin2-iBeg2
1448                        dMdpk2 = np.zeros(shape=(6,lenBF2))
1449                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
1450                        for i in range(5):
1451                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[13]*refl[9]*kRatio*dMdipk2[i]
1452                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[13]*dMdipk2[0]
1453                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1454                if Phase['General'].get('doPawley'):
1455                    dMdpw = np.zeros(len(x))
1456                    try:
1457                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1458                        idx = varylist.index(pIdx)
1459                        dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1460                        if Ka2:
1461                            dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1462                        dMdv[idx] = dMdpw
1463                    except: # ValueError:
1464                        pass
1465                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1466                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1467                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1468                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1469                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1470                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1471                    hfx+'DisplaceY':[dpdY,'pos'],hfx+'Absorption':[dFdAb,'int'],}
1472                for name in names:
1473                    item = names[name]
1474                    if name in varylist:
1475                        dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1476                        if Ka2:
1477                            dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1478                    elif name in dependentVars:
1479                        depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1480                        if Ka2:
1481                            depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1482                for iPO in dIdPO:
1483                    if iPO in varylist:
1484                        dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1485                        if Ka2:
1486                            dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1487                    elif iPO in dependentVars:
1488                        depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1489                        if Ka2:
1490                            depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1491                for i,name in enumerate(['omega','chi','phi']):
1492                    aname = pfx+'SH '+name
1493                    if aname in varylist:
1494                        dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1495                        if Ka2:
1496                            dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1497                    elif aname in dependentVars:
1498                        depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1499                        if Ka2:
1500                            depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1501                for iSH in dFdODF:
1502                    if iSH in varylist:
1503                        dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1504                        if Ka2:
1505                            dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1506                    elif iSH in dependentVars:
1507                        depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1508                        if Ka2:
1509                            depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1510                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1511                for name,dpdA in cellDervNames:
1512                    if name in varylist:
1513                        dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1514                        if Ka2:
1515                            dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1516                    elif name in dependentVars:
1517                        depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1518                        if Ka2:
1519                            depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1520                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1521                for name in dDijDict:
1522                    if name in varylist:
1523                        dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1524                        if Ka2:
1525                            dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1526                    elif name in dependentVars:
1527                        depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1528                        if Ka2:
1529                            depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1530                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
1531                for name in gamDict:
1532                    if name in varylist:
1533                        dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1534                        if Ka2:
1535                            dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1536                    elif name in dependentVars:
1537                        depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1538                        if Ka2:
1539                            depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1540                for name in sigDict:
1541                    if name in varylist:
1542                        dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1543                        if Ka2:
1544                            dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1545                    elif name in dependentVars:
1546                        depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1547                        if Ka2:
1548                            depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1549                for name in ['BabA','BabU']:
1550                    if phfx+name in varylist:
1551                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
1552                        if Ka2:
1553                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]
1554                    elif phfx+name in dependentVars:                   
1555                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
1556                        if Ka2:
1557                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]                 
1558            elif 'T' in calcControls[hfx+'histType']:
1559                print 'TOF Undefined at present'
1560                raise Exception    #no TOF yet
1561            if not Phase['General'].get('doPawley'):
1562                #do atom derivatives -  for RB,F,X & U so far             
1563                corr = dervDict['int']/refl[9]
1564                if Ka2:
1565                    corr2 = dervDict2['int']/refl[9]
1566                for name in varylist+dependentVars:
1567                    if '::RBV;' in name:
1568                        pass
1569                    else:
1570                        try:
1571                            aname = name.split(pfx)[1][:2]
1572                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1573                        except IndexError:
1574                            continue
1575                    if name in varylist:
1576                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1577                        if Ka2:
1578                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1579                    elif name in dependentVars:
1580                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1581                        if Ka2:
1582                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1583#        print 'profile derv time: %.3fs'%(time.time()-time0)
1584    # now process derivatives in constraints
1585    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1586    return dMdv
1587
1588def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1589    '''Loop over histograms and compute derivatives of the fitting
1590    model (M) with respect to all parameters.  Results are returned in
1591    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
1592    n is the number of parameters and m is the number of data
1593    points. This can exceed memory when m gets large. This routine is
1594    used when refinement derivatives are selected as "analtytic
1595    Jacobian" in Controls.
1596
1597    :returns: Jacobian numpy.array dMdv for all histograms concatinated
1598    '''
1599    parmDict.update(zip(varylist,values))
1600    G2mv.Dict2Map(parmDict,varylist)
1601    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1602    nvar = len(varylist)
1603    dMdv = np.empty(0)
1604    histoList = Histograms.keys()
1605    histoList.sort()
1606    for histogram in histoList:
1607        if 'PWDR' in histogram[:4]:
1608            Histogram = Histograms[histogram]
1609            hId = Histogram['hId']
1610            hfx = ':%d:'%(hId)
1611            wtFactor = calcControls[hfx+'wtFactor']
1612            Limits = calcControls[hfx+'Limits']
1613            x,y,w,yc,yb,yd = Histogram['Data']
1614            W = wtFactor*w
1615            xB = np.searchsorted(x,Limits[0])
1616            xF = np.searchsorted(x,Limits[1])
1617            dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
1618                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1619        elif 'HKLF' in histogram[:4]:
1620            Histogram = Histograms[histogram]
1621            nobs = Histogram['Nobs']
1622            phase = Histogram['Reflection Lists']
1623            Phase = Phases[phase]
1624            hId = Histogram['hId']
1625            hfx = ':%d:'%(hId)
1626            wtFactor = calcControls[hfx+'wtFactor']
1627            pfx = '%d::'%(Phase['pId'])
1628            phfx = '%d:%d:'%(Phase['pId'],hId)
1629            SGData = Phase['General']['SGData']
1630            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1631            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1632            refList = Histogram['Data']
1633            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1634            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1635            dMdvh = np.zeros((len(varylist),len(refList)))
1636            dependentVars = G2mv.GetDependentVars()
1637            depDerivDict = {}
1638            for j in dependentVars:
1639                depDerivDict[j] = np.zeros(shape=(len(refList)))
1640            if calcControls['F**2']:
1641                for iref,ref in enumerate(refList):
1642                    if ref[6] > 0:
1643                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1644                        w = 1.0/ref[6]
1645                        if w*ref[5] >= calcControls['minF/sig']:
1646                            for j,var in enumerate(varylist):
1647                                if var in dFdvDict:
1648                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor
1649                            for var in dependentVars:
1650                                if var in dFdvDict:
1651                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor
1652                            if phfx+'Scale' in varylist:
1653                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1654                            elif phfx+'Scale' in dependentVars:
1655                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor
1656                            for item in ['Ep','Es','Eg']:
1657                                if phfx+item in varylist:
1658                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1659                                elif phfx+item in dependentVars:
1660                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1661                            for item in ['BabA','BabU']:
1662                                if phfx+item in varylist:
1663                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1664                                elif phfx+item in dependentVars:
1665                                    depDerivDict[phfx+item][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1666            else:
1667                for iref,ref in enumerate(refList):
1668                    if ref[5] > 0.:
1669                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1670                        Fo = np.sqrt(ref[5])
1671                        Fc = np.sqrt(ref[7])
1672                        w = 1.0/ref[6]
1673                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1674                            for j,var in enumerate(varylist):
1675                                if var in dFdvDict:
1676                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1677                            for var in dependentVars:
1678                                if var in dFdvDict:
1679                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1680                            if phfx+'Scale' in varylist:
1681                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1682                            elif phfx+'Scale' in dependentVars:
1683                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor                           
1684                            for item in ['Ep','Es','Eg']:
1685                                if phfx+item in varylist:
1686                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1687                                elif phfx+item in dependentVars:
1688                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1689                            for item in ['BabA','BabU']:
1690                                if phfx+item in varylist:
1691                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1692                                elif phfx+item in dependentVars:
1693                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1694            # now process derivatives in constraints
1695            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1696        else:
1697            continue        #skip non-histogram entries
1698        if len(dMdv):
1699            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
1700        else:
1701            dMdv = np.sqrt(wtFactor)*dMdvh
1702           
1703    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1704    if np.any(pVals):
1705        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1706        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
1707       
1708    return dMdv
1709
1710def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1711    '''Loop over histograms and compute derivatives of the fitting
1712    model (M) with respect to all parameters.  For each histogram, the
1713    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
1714    number of parameters and m is the number of data points *in the
1715    histogram*. The (n by n) Hessian is computed from each Jacobian
1716    and it is returned.  This routine is used when refinement
1717    derivatives are selected as "analtytic Hessian" in Controls.
1718
1719    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
1720    '''
1721    'Needs a doc string'
1722    parmDict.update(zip(varylist,values))
1723    G2mv.Dict2Map(parmDict,varylist)
1724    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1725    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
1726    nvar = len(varylist)
1727    Hess = np.empty(0)
1728    histoList = Histograms.keys()
1729    histoList.sort()
1730    for histogram in histoList:
1731        if 'PWDR' in histogram[:4]:
1732            Histogram = Histograms[histogram]
1733            hId = Histogram['hId']
1734            hfx = ':%d:'%(hId)
1735            wtFactor = calcControls[hfx+'wtFactor']
1736            Limits = calcControls[hfx+'Limits']
1737            x,y,w,yc,yb,yd = Histogram['Data']
1738            W = wtFactor*w
1739            dy = y-yc
1740            xB = np.searchsorted(x,Limits[0])
1741            xF = np.searchsorted(x,Limits[1])
1742            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
1743                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1744            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
1745            Dy = dy[xB:xF][np.newaxis,:]
1746            dMdvh *= Wt
1747            if dlg:
1748                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1749            if len(Hess):
1750                Hess += np.inner(dMdvh,dMdvh)
1751                dMdvh *= Wt*Dy
1752                Vec += np.sum(dMdvh,axis=1)
1753            else:
1754                Hess = np.inner(dMdvh,dMdvh)
1755                dMdvh *= Wt*Dy
1756                Vec = np.sum(dMdvh,axis=1)
1757        elif 'HKLF' in histogram[:4]:
1758            Histogram = Histograms[histogram]
1759            nobs = Histogram['Nobs']
1760            phase = Histogram['Reflection Lists']
1761            Phase = Phases[phase]
1762            hId = Histogram['hId']
1763            hfx = ':%d:'%(hId)
1764            wtFactor = calcControls[hfx+'wtFactor']
1765            pfx = '%d::'%(Phase['pId'])
1766            phfx = '%d:%d:'%(Phase['pId'],hId)
1767            SGData = Phase['General']['SGData']
1768            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1769            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1770            refList = Histogram['Data']
1771            time0 = time.time()
1772            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)  #accurate for powders!
1773            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1774            dMdvh = np.zeros((len(varylist),len(refList)))
1775            dependentVars = G2mv.GetDependentVars()
1776            depDerivDict = {}
1777            for j in dependentVars:
1778                depDerivDict[j] = np.zeros(shape=(len(refList)))
1779            wdf = np.zeros(len(refList))
1780            if calcControls['F**2']:
1781                for iref,ref in enumerate(refList):
1782                    if ref[6] > 0:
1783                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1784                        w =  1.0/ref[6]
1785                        if w*ref[5] >= calcControls['minF/sig']:
1786                            wdf[iref] = w*(ref[5]-ref[7])
1787                            for j,var in enumerate(varylist):
1788                                if var in dFdvDict:
1789                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1790                            for var in dependentVars:
1791                                if var in dFdvDict:
1792                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1793                            if phfx+'Scale' in varylist:
1794                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1795                            elif phfx+'Scale' in dependentVars:
1796                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor
1797                            for item in ['Ep','Es','Eg']:
1798                                if phfx+item in varylist:
1799                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1800                                elif phfx+item in dependentVars:
1801                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1802                            for item in ['BabA','BabU']:
1803                                if phfx+item in varylist:
1804                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1805                                elif phfx+item in dependentVars:
1806                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1807            else:
1808                for iref,ref in enumerate(refList):
1809                    if ref[5] > 0.:
1810                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1811                        Fo = np.sqrt(ref[5])
1812                        Fc = np.sqrt(ref[7])
1813                        w = 1.0/ref[6]
1814                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1815                            wdf[iref] = 2.0*Fo*w*(Fo-Fc)
1816                            for j,var in enumerate(varylist):
1817                                if var in dFdvDict:
1818                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1819                            for var in dependentVars:
1820                                if var in dFdvDict:
1821                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1822                            if phfx+'Scale' in varylist:
1823                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1824                            elif phfx+'Scale' in dependentVars:
1825                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor                           
1826                            for item in ['Ep','Es','Eg']:
1827                                if phfx+item in varylist:
1828                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1829                                elif phfx+item in dependentVars:
1830                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1831                            for item in ['BabA','BabU']:
1832                                if phfx+item in varylist:
1833                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1834                                elif phfx+item in dependentVars:
1835                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1836            # now process derivatives in constraints
1837            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1838
1839            if dlg:
1840                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1841            if len(Hess):
1842                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
1843                Hess += wtFactor*np.inner(dMdvh,dMdvh)
1844            else:
1845                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
1846                Hess = wtFactor*np.inner(dMdvh,dMdvh)
1847        else:
1848            continue        #skip non-histogram entries
1849    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1850    if np.any(pVals):
1851        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1852        Vec += np.sum(dpdv*pWt*pVals,axis=1)
1853        Hess += np.inner(dpdv*pWt,dpdv)
1854    return Vec,Hess
1855
1856def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):       
1857    'Needs a doc string'
1858    parmDict.update(zip(varylist,values))
1859    Values2Dict(parmDict, varylist, values)
1860    G2mv.Dict2Map(parmDict,varylist)
1861    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1862    M = np.empty(0)
1863    SumwYo = 0
1864    Nobs = 0
1865    ApplyRBModels(parmDict,Phases,rigidbodyDict)
1866    histoList = Histograms.keys()
1867    histoList.sort()
1868    for histogram in histoList:
1869        if 'PWDR' in histogram[:4]:
1870            Histogram = Histograms[histogram]
1871            hId = Histogram['hId']
1872            hfx = ':%d:'%(hId)
1873            wtFactor = calcControls[hfx+'wtFactor']
1874            Limits = calcControls[hfx+'Limits']
1875            x,y,w,yc,yb,yd = Histogram['Data']
1876            W = wtFactor*w
1877            yc *= 0.0                           #zero full calcd profiles
1878            yb *= 0.0
1879            yd *= 0.0
1880            xB = np.searchsorted(x,Limits[0])
1881            xF = np.searchsorted(x,Limits[1])
1882            Histogram['Nobs'] = ma.count(x[xB:xF])
1883            Nobs += Histogram['Nobs']
1884            Histogram['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
1885            SumwYo += Histogram['sumwYo']
1886            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
1887                varylist,Histogram,Phases,calcControls,pawleyLookup)
1888            yc[xB:xF] += yb[xB:xF]
1889            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
1890            Histogram['sumwYd'] = ma.sum(np.sqrt(W[xB:xF])*(yd[xB:xF]))
1891            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
1892            Histogram['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['sumwYo'])*100.)
1893            if dlg:
1894                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1895            M = np.concatenate((M,wdy))
1896#end of PWDR processing
1897        elif 'HKLF' in histogram[:4]:
1898            Histogram = Histograms[histogram]
1899            phase = Histogram['Reflection Lists']
1900            Phase = Phases[phase]
1901            hId = Histogram['hId']
1902            hfx = ':%d:'%(hId)
1903            wtFactor = calcControls[hfx+'wtFactor']
1904            pfx = '%d::'%(Phase['pId'])
1905            phfx = '%d:%d:'%(Phase['pId'],hId)
1906            SGData = Phase['General']['SGData']
1907            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1908            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1909            refList = Histogram['Data']
1910            time0 = time.time()
1911            StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1912#            print 'sf-calc time: %.3f'%(time.time()-time0)
1913            df = np.zeros(len(refList))
1914            sumwYo = 0
1915            sumFo = 0
1916            sumFo2 = 0
1917            sumdF = 0
1918            sumdF2 = 0
1919            nobs = 0
1920            if calcControls['F**2']:
1921                for i,ref in enumerate(refList):
1922                    if ref[6] > 0:
1923                        SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1924                        w = 1.0/ref[6]
1925                        ref[7] = parmDict[phfx+'Scale']*ref[9]
1926                        ref[7] *= ref[13]                       #correct Fc^2 for extinction
1927                        ref[8] = ref[5]/parmDict[phfx+'Scale']
1928                        if w*ref[5] >= calcControls['minF/sig']:
1929                            sumFo2 += ref[5]
1930                            Fo = np.sqrt(ref[5])
1931                            sumFo += Fo
1932                            sumFo2 += ref[5]
1933                            sumdF += abs(Fo-np.sqrt(ref[7]))
1934                            sumdF2 += abs(ref[5]-ref[7])
1935                            nobs += 1
1936                            df[i] = -w*(ref[5]-ref[7])
1937                            sumwYo += (w*ref[5])**2
1938            else:
1939                for i,ref in enumerate(refList):
1940                    if ref[5] > 0.:
1941                        SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1942                        ref[7] = parmDict[phfx+'Scale']*ref[9]
1943                        ref[7] *= ref[13]                       #correct Fc^2 for extinction
1944                        Fo = np.sqrt(ref[5])
1945                        Fc = np.sqrt(ref[7])
1946                        w = 2.0*Fo/ref[6]
1947                        if w*Fo >= calcControls['minF/sig']:
1948                            sumFo += Fo
1949                            sumFo2 += ref[5]
1950                            sumdF += abs(Fo-Fc)
1951                            sumdF2 += abs(ref[5]-ref[7])
1952                            nobs += 1
1953                            df[i] = -w*(Fo-Fc)
1954                            sumwYo += (w*Fo)**2
1955            Histogram['Nobs'] = nobs
1956            Histogram['sumwYo'] = sumwYo
1957            SumwYo += sumwYo
1958            Histogram['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['sumwYo'])*100.)
1959            Histogram[phfx+'Rf'] = 100.*sumdF/sumFo
1960            Histogram[phfx+'Rf^2'] = 100.*sumdF2/sumFo2
1961            Histogram[phfx+'Nref'] = nobs
1962            Nobs += nobs
1963            if dlg:
1964                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1965            M = np.concatenate((M,wtFactor*df))
1966# end of HKLF processing
1967    Histograms['sumwYo'] = SumwYo
1968    Histograms['Nobs'] = Nobs
1969    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
1970    if dlg:
1971        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
1972        if not GoOn:
1973            parmDict['saved values'] = values
1974            dlg.Destroy()
1975            raise Exception         #Abort!!
1976    pDict,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1977    if np.any(pVals):
1978        pSum = np.sum(pWt*pVals**2)
1979        print 'Penalty function: %.3f on %d terms'%(pSum,len(pVals))
1980        Nobs += len(pVals)
1981        M = np.concatenate((M,np.sqrt(pWt)*pVals))
1982    return M
1983                       
Note: See TracBrowser for help on using the repository browser.