source: branch/2frame/GSASIIstrMath.py @ 2947

Last change on this file since 2947 was 2947, checked in by toby, 4 years ago

check in mp fix before reverting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 201.7 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2017-07-30 00:10:21 +0000 (Sun, 30 Jul 2017) $
8# $Author: toby $
9# $Revision: 2947 $
10# $URL: branch/2frame/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 2947 2017-07-30 00:10:21Z toby $
12########### SVN repository information ###################
13import time
14import copy
15import multiprocessing as mp
16import numpy as np
17import numpy.ma as ma
18import numpy.linalg as nl
19import scipy.stats as st
20import GSASIIpath
21GSASIIpath.SetVersionNumber("$Revision: 2947 $")
22import GSASIIElem as G2el
23import GSASIIlattice as G2lat
24import GSASIIspc as G2spc
25import GSASIIpwd as G2pwd
26import GSASIImapvars as G2mv
27import GSASIImath as G2mth
28import GSASIIobj as G2obj
29import GSASIImpsubs as G2mp
30
31sind = lambda x: np.sin(x*np.pi/180.)
32cosd = lambda x: np.cos(x*np.pi/180.)
33tand = lambda x: np.tan(x*np.pi/180.)
34asind = lambda x: 180.*np.arcsin(x)/np.pi
35acosd = lambda x: 180.*np.arccos(x)/np.pi
36atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
37   
38ateln2 = 8.0*np.log(2.0)
39twopi = 2.0*np.pi
40twopisq = 2.0*np.pi**2
41nxs = np.newaxis
42
43################################################################################
44##### Rigid Body Models
45################################################################################
46       
47def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
48    ''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
49    current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
50    '''
51    atxIds = ['Ax:','Ay:','Az:']
52    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
53    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
54    if not RBIds['Vector'] and not RBIds['Residue']:
55        return
56    VRBIds = RBIds['Vector']
57    RRBIds = RBIds['Residue']
58    if Update:
59        RBData = rigidbodyDict
60    else:
61        RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
62    if RBIds['Vector']:                       # first update the vector magnitudes
63        VRBData = RBData['Vector']
64        for i,rbId in enumerate(VRBIds):
65            if VRBData[rbId]['useCount']:
66                for j in range(len(VRBData[rbId]['VectMag'])):
67                    name = '::RBV;'+str(j)+':'+str(i)
68                    VRBData[rbId]['VectMag'][j] = parmDict[name]
69    for phase in Phases:
70        Phase = Phases[phase]
71        General = Phase['General']
72        cx,ct,cs,cia = General['AtomPtrs']
73        cell = General['Cell'][1:7]
74        Amat,Bmat = G2lat.cell2AB(cell)
75        AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
76        pfx = str(Phase['pId'])+'::'
77        if Update:
78            RBModels = Phase['RBModels']
79        else:
80            RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
81        for irb,RBObj in enumerate(RBModels.get('Vector',[])):
82            jrb = VRBIds.index(RBObj['RBId'])
83            rbsx = str(irb)+':'+str(jrb)
84            for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
85                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
86            for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
87                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
88            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
89            TLS = RBObj['ThermalMotion']
90            if 'T' in TLS[0]:
91                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
92                    TLS[1][i] = parmDict[pfx+pt+rbsx]
93            if 'L' in TLS[0]:
94                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
95                    TLS[1][i+6] = parmDict[pfx+pt+rbsx]
96            if 'S' in TLS[0]:
97                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
98                    TLS[1][i+12] = parmDict[pfx+pt+rbsx]
99            if 'U' in TLS[0]:
100                TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
101            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
102            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
103            for i,x in enumerate(XYZ):
104                atId = RBObj['Ids'][i]
105                for j in [0,1,2]:
106                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
107                if UIJ[i][0] == 'A':
108                    for j in range(6):
109                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
110                elif UIJ[i][0] == 'I':
111                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
112           
113        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
114            jrb = RRBIds.index(RBObj['RBId'])
115            rbsx = str(irb)+':'+str(jrb)
116            for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
117                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
118            for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
119                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]               
120            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
121            TLS = RBObj['ThermalMotion']
122            if 'T' in TLS[0]:
123                for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
124                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
125            if 'L' in TLS[0]:
126                for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
127                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
128            if 'S' in TLS[0]:
129                for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
130                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
131            if 'U' in TLS[0]:
132                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
133            for itors,tors in enumerate(RBObj['Torsions']):
134                tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
135            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
136            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
137            for i,x in enumerate(XYZ):
138                atId = RBObj['Ids'][i]
139                for j in [0,1,2]:
140                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
141                if UIJ[i][0] == 'A':
142                    for j in range(6):
143                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
144                elif UIJ[i][0] == 'I':
145                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
146                   
147def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
148    'Needs a doc string'
149    atxIds = ['dAx:','dAy:','dAz:']
150    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
151    OIds = ['Oa:','Oi:','Oj:','Ok:']
152    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
153    if not RBIds['Vector'] and not RBIds['Residue']:
154        return
155    VRBIds = RBIds['Vector']
156    RRBIds = RBIds['Residue']
157    RBData = rigidbodyDict
158    for item in parmDict:
159        if 'RB' in item:
160            dFdvDict[item] = 0.        #NB: this is a vector which is no. refl. long & must be filled!
161    General = Phase['General']
162    cx,ct,cs,cia = General['AtomPtrs']
163    cell = General['Cell'][1:7]
164    Amat,Bmat = G2lat.cell2AB(cell)
165    rpd = np.pi/180.
166    rpd2 = rpd**2
167    g = nl.inv(np.inner(Bmat,Bmat))
168    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
169        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
170    AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
171    pfx = str(Phase['pId'])+'::'
172    RBModels =  Phase['RBModels']
173    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
174        VModel = RBData['Vector'][RBObj['RBId']]
175        Q = RBObj['Orient'][0]
176        jrb = VRBIds.index(RBObj['RBId'])
177        rbsx = str(irb)+':'+str(jrb)
178        dXdv = []
179        for iv in range(len(VModel['VectMag'])):
180            dCdv = []
181            for vec in VModel['rbVect'][iv]:
182                dCdv.append(G2mth.prodQVQ(Q,vec))
183            dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
184        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
185        for ia,atId in enumerate(RBObj['Ids']):
186            atNum = AtLookup[atId]
187            dx = 0.00001
188            for iv in range(len(VModel['VectMag'])):
189                for ix in [0,1,2]:
190                    dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
191            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
192                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
193            for iv in range(4):
194                Q[iv] -= dx
195                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
196                Q[iv] += 2.*dx
197                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
198                Q[iv] -= dx
199                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
200                for ix in [0,1,2]:
201                    dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
202            X = G2mth.prodQVQ(Q,Cart[ia])
203            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
204            dFdu = G2lat.U6toUij(dFdu.T)
205            dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))           
206            dFdu = G2lat.UijtoU6(dFdu)
207            atNum = AtLookup[atId]
208            if 'T' in RBObj['ThermalMotion'][0]:
209                for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
210                    dFdvDict[pfx+name+rbsx] += dFdu[i]
211            if 'L' in RBObj['ThermalMotion'][0]:
212                dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
213                dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
214                dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
215                dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
216                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
217                dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2*(-dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
218                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
219                dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2*(-dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
220                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
221            if 'S' in RBObj['ThermalMotion'][0]:
222                dFdvDict[pfx+'RBVS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
223                dFdvDict[pfx+'RBVS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
224                dFdvDict[pfx+'RBVS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
225                dFdvDict[pfx+'RBVS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
226                dFdvDict[pfx+'RBVS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
227                dFdvDict[pfx+'RBVS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
228                dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
229                dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
230            if 'U' in RBObj['ThermalMotion'][0]:
231                dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
232
233
234    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
235        Q = RBObj['Orient'][0]
236        jrb = RRBIds.index(RBObj['RBId'])
237        torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
238        rbsx = str(irb)+':'+str(jrb)
239        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
240        for itors,tors in enumerate(RBObj['Torsions']):     #derivative error?
241            tname = pfx+'RBRTr;'+str(itors)+':'+rbsx           
242            orId,pvId = torData[itors][:2]
243            pivotVec = Cart[orId]-Cart[pvId]
244            QA = G2mth.AVdeg2Q(-0.001,pivotVec)
245            QB = G2mth.AVdeg2Q(0.001,pivotVec)
246            for ir in torData[itors][3]:
247                atNum = AtLookup[RBObj['Ids'][ir]]
248                rVec = Cart[ir]-Cart[pvId]
249                dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
250                dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
251                for ix in [0,1,2]:
252                    dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
253        for ia,atId in enumerate(RBObj['Ids']):
254            atNum = AtLookup[atId]
255            dx = 0.00001
256            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
257                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
258            for iv in range(4):
259                Q[iv] -= dx
260                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
261                Q[iv] += 2.*dx
262                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
263                Q[iv] -= dx
264                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
265                for ix in [0,1,2]:
266                    dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
267            X = G2mth.prodQVQ(Q,Cart[ia])
268            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
269            dFdu = G2lat.U6toUij(dFdu.T)
270            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
271            dFdu = G2lat.UijtoU6(dFdu)
272            atNum = AtLookup[atId]
273            if 'T' in RBObj['ThermalMotion'][0]:
274                for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
275                    dFdvDict[pfx+name+rbsx] += dFdu[i]
276            if 'L' in RBObj['ThermalMotion'][0]:
277                dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
278                dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
279                dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
280                dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
281                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
282                dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2*(dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
283                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
284                dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2*(dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
285                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
286            if 'S' in RBObj['ThermalMotion'][0]:
287                dFdvDict[pfx+'RBRS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
288                dFdvDict[pfx+'RBRS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
289                dFdvDict[pfx+'RBRS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
290                dFdvDict[pfx+'RBRS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
291                dFdvDict[pfx+'RBRS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
292                dFdvDict[pfx+'RBRS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
293                dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
294                dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
295            if 'U' in RBObj['ThermalMotion'][0]:
296                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
297   
298################################################################################
299##### Penalty & restraint functions
300################################################################################
301
302def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
303    'Needs a doc string'
304    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
305    pNames = []
306    pVals = []
307    pWt = []
308    negWt = {}
309    pWsum = {}
310    pWnum = {}
311    for phase in Phases:
312        pId = Phases[phase]['pId']
313        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
314        General = Phases[phase]['General']
315        cx,ct,cs,cia = General['AtomPtrs']
316        textureData = General['SH Texture']
317        SGData = General['SGData']
318        Atoms = Phases[phase]['Atoms']
319        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
320        cell = General['Cell'][1:7]
321        Amat,Bmat = G2lat.cell2AB(cell)
322        if phase not in restraintDict:
323            continue
324        phaseRest = restraintDict[phase]
325        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
326            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
327            ['ChemComp','Sites'],['Texture','HKLs'],]
328        for name,rest in names:
329            pWsum[name] = 0.
330            pWnum[name] = 0
331            itemRest = phaseRest[name]
332            if itemRest[rest] and itemRest['Use']:
333                wt = itemRest['wtFactor']
334                if name in ['Bond','Angle','Plane','Chiral']:
335                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
336                        pNames.append(str(pId)+':'+name+':'+str(i))
337                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
338                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
339                        if name == 'Bond':
340                            calc = G2mth.getRestDist(XYZ,Amat)
341                        elif name == 'Angle':
342                            calc = G2mth.getRestAngle(XYZ,Amat)
343                        elif name == 'Plane':
344                            calc = G2mth.getRestPlane(XYZ,Amat)
345                        elif name == 'Chiral':
346                            calc = G2mth.getRestChiral(XYZ,Amat)
347                        pVals.append(obs-calc)
348                        pWt.append(wt/esd**2)
349                        pWsum[name] += wt*((obs-calc)/esd)**2
350                        pWnum[name] += 1
351                elif name in ['Torsion','Rama']:
352                    coeffDict = itemRest['Coeff']
353                    for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
354                        pNames.append(str(pId)+':'+name+':'+str(i))
355                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
356                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
357                        if name == 'Torsion':
358                            tor = G2mth.getRestTorsion(XYZ,Amat)
359                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
360                        else:
361                            phi,psi = G2mth.getRestRama(XYZ,Amat)
362                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
363                        pVals.append(restr)
364                        pWt.append(wt/esd**2)
365                        pWsum[name] += wt*(restr/esd)**2
366                        pWnum[name] += 1
367                elif name == 'ChemComp':
368                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
369                        pNames.append(str(pId)+':'+name+':'+str(i))
370                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
371                        frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
372                        calc = np.sum(mul*frac*factors)
373                        pVals.append(obs-calc)
374                        pWt.append(wt/esd**2)                   
375                        pWsum[name] += wt*((obs-calc)/esd)**2
376                        pWnum[name] += 1
377                elif name == 'Texture':
378                    SHkeys = textureData['SH Coeff'][1].keys()
379                    SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
380                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
381                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
382                    for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
383                        PH = np.array(hkl)
384                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
385                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
386                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
387                        Z1 = ma.masked_greater(Z,0.0)           #is this + or -?
388                        IndZ1 = np.array(ma.nonzero(Z1))
389                        for ind in IndZ1.T:
390                            pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
391                            pVals.append(Z1[ind[0]][ind[1]])
392                            pWt.append(wt/esd1**2)
393                            pWsum[name] += wt*(-Z1[ind[0]][ind[1]]/esd1)**2
394                            pWnum[name] += 1
395                        if ifesd2:
396                            Z2 = 1.-Z
397                            for ind in np.ndindex(grid,grid):
398                                pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
399                                pVals.append(Z2[ind[0]][ind[1]])
400                                pWt.append(wt/esd2**2)
401                                pWsum[name] += wt*(Z2/esd2)**2
402                                pWnum[name] += 1
403       
404    for phase in Phases:
405        name = 'SH-Pref.Ori.'
406        pId = Phases[phase]['pId']
407        General = Phases[phase]['General']
408        SGData = General['SGData']
409        cell = General['Cell'][1:7]
410        pWsum[name] = 0.0
411        pWnum[name] = 0
412        for hist in Phases[phase]['Histograms']:
413            if not Phases[phase]['Histograms'][hist]['Use']:
414                continue
415            if hist in Histograms and 'PWDR' in hist:
416                hId = Histograms[hist]['hId']
417                phfx = '%d:%d:'%(pId,hId)
418                if calcControls[phfx+'poType'] == 'SH':
419                    toler = calcControls[phfx+'SHtoler']
420                    wt = 1./toler**2
421                    HKLs = np.array(calcControls[phfx+'SHhkl'])
422                    SHnames = calcControls[phfx+'SHnames']
423                    SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
424                    for i,PH in enumerate(HKLs):
425                        phi,beta = G2lat.CrsAng(PH,cell,SGData)
426                        SH3Coef = {}
427                        for item in SHcof:
428                            L,N = eval(item.strip('C'))
429                            SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
430                        ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
431                        X = np.linspace(0,90.0,26)
432                        Y = ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0)       #+ or -?
433                        IndY = ma.nonzero(Y)
434                        for ind in IndY[0]:
435                            pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind]))
436                            pVals.append(Y[ind])
437                            pWt.append(wt)
438                            pWsum[name] += wt*(Y[ind])**2
439                            pWnum[name] += 1
440    pWsum['PWLref'] = 0.
441    pWnum['PWLref'] = 0
442    for item in varyList:
443        if 'PWLref' in item and parmDict[item] < 0.:
444            pId = int(item.split(':')[0])
445            if negWt[pId]:
446                pNames.append(item)
447                pVals.append(parmDict[item])
448                pWt.append(negWt[pId])
449                pWsum['PWLref'] += negWt[pId]*(parmDict[item])**2
450                pWnum['PWLref'] += 1
451    pVals = np.array(pVals)
452    pWt = np.array(pWt)         #should this be np.sqrt?
453    return pNames,pVals,pWt,pWsum,pWnum
454   
455def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
456    'Needs a doc string'
457    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
458    pDerv = np.zeros((len(varyList),len(pVal)))
459    for phase in Phases:
460#        if phase not in restraintDict:
461#            continue
462        pId = Phases[phase]['pId']
463        General = Phases[phase]['General']
464        cx,ct,cs,cia = General['AtomPtrs']
465        SGData = General['SGData']
466        Atoms = Phases[phase]['Atoms']
467        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
468        cell = General['Cell'][1:7]
469        Amat,Bmat = G2lat.cell2AB(cell)
470        textureData = General['SH Texture']
471
472        SHkeys = textureData['SH Coeff'][1].keys()
473        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
474        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
475        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
476        sam = SamSym[textureData['Model']]
477        phaseRest = restraintDict.get(phase,{})
478        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
479            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
480            'ChemComp':'Sites','Texture':'HKLs'}
481        lasthkl = np.array([0,0,0])
482        for ip,pName in enumerate(pNames):
483            pnames = pName.split(':')
484            if pId == int(pnames[0]):
485                name = pnames[1]
486                if 'PWL' in pName:
487                    pDerv[varyList.index(pName)][ip] += 1.
488                    continue
489                elif 'SH-' in pName:
490                    continue
491                id = int(pnames[2]) 
492                itemRest = phaseRest[name]
493                if name in ['Bond','Angle','Plane','Chiral']:
494                    indx,ops,obs,esd = itemRest[names[name]][id]
495                    dNames = []
496                    for ind in indx:
497                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
498                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
499                    if name == 'Bond':
500                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
501                    elif name == 'Angle':
502                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
503                    elif name == 'Plane':
504                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
505                    elif name == 'Chiral':
506                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
507                elif name in ['Torsion','Rama']:
508                    coffDict = itemRest['Coeff']
509                    indx,ops,cofName,esd = itemRest[names[name]][id]
510                    dNames = []
511                    for ind in indx:
512                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
513                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
514                    if name == 'Torsion':
515                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
516                    else:
517                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
518                elif name == 'ChemComp':
519                    indx,factors,obs,esd = itemRest[names[name]][id]
520                    dNames = []
521                    for ind in indx:
522                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
523                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
524                        deriv = mul*factors
525                elif 'Texture' in name:
526                    deriv = []
527                    dNames = []
528                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
529                    hkl = np.array(hkl)
530                    if np.any(lasthkl-hkl):
531                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
532                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
533                        lasthkl = copy.copy(hkl)                       
534                    if 'unit' in name:
535                        pass
536                    else:
537                        gam = float(pnames[3])
538                        psi = float(pnames[4])
539                        for SHname in ODFln:
540                            l,m,n = eval(SHname[1:])
541                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
542                            dNames += [str(pId)+'::'+SHname]
543                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
544                for dName,drv in zip(dNames,deriv):
545                    try:
546                        ind = varyList.index(dName)
547                        pDerv[ind][ip] += drv
548                    except ValueError:
549                        pass
550       
551        lasthkl = np.array([0,0,0])
552        for ip,pName in enumerate(pNames):
553            deriv = []
554            dNames = []
555            pnames = pName.split(':')
556            if 'SH-' in pName and pId == int(pnames[0]):
557                hId = int(pnames[1])
558                phfx = '%d:%d:'%(pId,hId)
559                psi = float(pnames[4])
560                HKLs = calcControls[phfx+'SHhkl']
561                SHnames = calcControls[phfx+'SHnames']
562                SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
563                hkl = np.array(HKLs[int(pnames[3])])     
564                if np.any(lasthkl-hkl):
565                    phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
566                    SH3Coef = {}
567                    for item in SHcof:
568                        L,N = eval(item.strip('C'))
569                        SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
570                    ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
571                    lasthkl = copy.copy(hkl)                       
572                for SHname in SHnames:
573                    l,n = eval(SHname[1:])
574                    SH3name = 'C%d,0,%d'%(l,n)
575                    Ksl = G2lat.GetKsl(l,0,'0',psi,0.0)[0]
576                    dNames += [phfx+SHname]
577                    deriv.append(ODFln[SH3name][0]*Ksl/SHcof[SHname])
578            for dName,drv in zip(dNames,deriv):
579                try:
580                    ind = varyList.index(dName)
581                    pDerv[ind][ip] += drv
582                except ValueError:
583                    pass
584    return pDerv
585
586################################################################################
587##### Function & derivative calculations
588################################################################################       
589                   
590def GetAtomFXU(pfx,calcControls,parmDict):
591    'Needs a doc string'
592    Natoms = calcControls['Natoms'][pfx]
593    Tdata = Natoms*[' ',]
594    Mdata = np.zeros(Natoms)
595    IAdata = Natoms*[' ',]
596    Fdata = np.zeros(Natoms)
597    Xdata = np.zeros((3,Natoms))
598    dXdata = np.zeros((3,Natoms))
599    Uisodata = np.zeros(Natoms)
600    Uijdata = np.zeros((6,Natoms))
601    Gdata = np.zeros((3,Natoms))
602    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
603        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
604        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
605        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
606        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5],
607        'AMx:':Gdata[0],'AMy:':Gdata[1],'AMz:':Gdata[2],}
608    for iatm in range(Natoms):
609        for key in keys:
610            parm = pfx+key+str(iatm)
611            if parm in parmDict:
612                keys[key][iatm] = parmDict[parm]
613    Fdata = np.where(Fdata,Fdata,1.e-8)         #avoid divide by zero in derivative calc.
614    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata
615   
616def GetAtomSSFXU(pfx,calcControls,parmDict):
617    'Needs a doc string'
618    Natoms = calcControls['Natoms'][pfx]
619    maxSSwave = calcControls['maxSSwave'][pfx]
620    Nwave = {'F':maxSSwave['Sfrac'],'X':maxSSwave['Spos'],'Y':maxSSwave['Spos'],'Z':maxSSwave['Spos'],
621        'U':maxSSwave['Sadp'],'M':maxSSwave['Smag'],'T':maxSSwave['Spos']}
622    XSSdata = np.zeros((6,maxSSwave['Spos'],Natoms))
623    FSSdata = np.zeros((2,maxSSwave['Sfrac'],Natoms))
624    USSdata = np.zeros((12,maxSSwave['Sadp'],Natoms))
625    MSSdata = np.zeros((6,maxSSwave['Smag'],Natoms))
626    waveTypes = []
627    keys = {'Fsin:':FSSdata[0],'Fcos:':FSSdata[1],'Fzero:':FSSdata[0],'Fwid:':FSSdata[1],
628        'Tmin:':XSSdata[0],'Tmax:':XSSdata[1],'Xmax:':XSSdata[2],'Ymax:':XSSdata[3],'Zmax:':XSSdata[4],
629        'Xsin:':XSSdata[0],'Ysin:':XSSdata[1],'Zsin:':XSSdata[2],'Xcos:':XSSdata[3],'Ycos:':XSSdata[4],'Zcos:':XSSdata[5],
630        'U11sin:':USSdata[0],'U22sin:':USSdata[1],'U33sin:':USSdata[2],'U12sin:':USSdata[3],'U13sin:':USSdata[4],'U23sin:':USSdata[5],
631        'U11cos:':USSdata[6],'U22cos:':USSdata[7],'U33cos:':USSdata[8],'U12cos:':USSdata[9],'U13cos:':USSdata[10],'U23cos:':USSdata[11],
632        'MXsin:':MSSdata[0],'MYsin:':MSSdata[1],'MZsin:':MSSdata[2],'MXcos:':MSSdata[3],'MYcos:':MSSdata[4],'MZcos:':MSSdata[5]}
633    for iatm in range(Natoms):
634        waveTypes.append(parmDict[pfx+'waveType:'+str(iatm)])
635        for key in keys:
636            for m in range(Nwave[key[0]]):
637                parm = pfx+key+str(iatm)+':%d'%(m)
638                if parm in parmDict:
639                    keys[key][m][iatm] = parmDict[parm]
640    return np.array(waveTypes),FSSdata,XSSdata,USSdata,MSSdata
641   
642def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
643    ''' Compute structure factors for all h,k,l for phase
644    puts the result, F^2, in each ref[8] in refList
645    operates on blocks of 100 reflections for speed
646    input:
647   
648    :param dict refDict: where
649        'RefList' list where each ref = h,k,l,it,d,...
650        'FF' dict of form factors - filed in below
651    :param np.array G:      reciprocal metric tensor
652    :param str pfx:    phase id string
653    :param dict SGData: space group info. dictionary output from SpcGroup
654    :param dict calcControls:
655    :param dict ParmDict:
656
657    '''       
658    phfx = pfx.split(':')[0]+hfx
659    ast = np.sqrt(np.diag(G))
660    Mast = twopisq*np.multiply.outer(ast,ast)
661    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
662    SGT = np.array([ops[1] for ops in SGData['SGOps']])
663    Ncen = len(SGData['SGCen'])
664    Nops = len(SGMT)*Ncen*(1+SGData['SGInv'])
665    FFtables = calcControls['FFtables']
666    BLtables = calcControls['BLtables']
667    MFtables = calcControls['MFtables']
668    Amat,Bmat = G2lat.Gmat2AB(G)
669    Flack = 1.0
670    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
671        Flack = 1.-2.*parmDict[phfx+'Flack']
672    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
673    TwDict = refDict.get('TwDict',{})           
674    if 'S' in calcControls[hfx+'histType']:
675        NTL = calcControls[phfx+'NTL']
676        NM = calcControls[phfx+'TwinNMN']+1
677        TwinLaw = calcControls[phfx+'TwinLaw']
678        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
679        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
680    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
681        GetAtomFXU(pfx,calcControls,parmDict)
682    if not Xdata.size:          #no atoms in phase!
683        return
684    if parmDict[pfx+'isMag']:
685        Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
686        Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
687        Gdata = np.inner(Bmat,Gdata.T)              #convert to crystal space
688        Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
689        if SGData['SGInv']:
690            Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
691        Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
692#        GSASIIpath.IPyBreak()
693        Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
694        Gdata = np.inner(Amat,Gdata.T)              #convert back to cart. space MXYZ, Natoms, NOps*Inv*Ncen
695        Gdata = np.swapaxes(Gdata,1,2)              # put Natoms last
696        Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T
697        if SGData['SGInv']:
698            Mag = np.repeat(Mag,2,axis=0)                  #Mag same shape as Gdata
699    if 'NC' in calcControls[hfx+'histType']:
700        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
701    elif 'X' in calcControls[hfx+'histType']:
702        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
703        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
704    Uij = np.array(G2lat.U6toUij(Uijdata))
705    bij = Mast*Uij.T
706    blkSize = 100       #no. of reflections in a block - size seems optimal
707    nRef = refDict['RefList'].shape[0]
708    SQ = 1./(2.*refDict['RefList'].T[4])**2
709    if 'N' in calcControls[hfx+'histType']:
710        dat = G2el.getBLvalues(BLtables)
711        refDict['FF']['El'] = dat.keys()
712        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
713        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
714        for iel,El in enumerate(refDict['FF']['El']):
715            if El in MFtables:
716                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
717    else:       #'X'
718        dat = G2el.getFFvalues(FFtables,0.)
719        refDict['FF']['El'] = dat.keys()
720        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
721        for iel,El in enumerate(refDict['FF']['El']):
722            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
723#reflection processing begins here - big arrays!
724    iBeg = 0
725#    time0 = time.time()
726    while iBeg < nRef:
727        iFin = min(iBeg+blkSize,nRef)
728        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
729        H = refl.T[:3]                          #array(blkSize,3)
730        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
731        TwMask = np.any(H,axis=-1)
732        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
733            for ir in range(blkSize):
734                iref = ir+iBeg
735                if iref in TwDict:
736                    for i in TwDict[iref]:
737                        for n in range(NTL):
738                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
739            TwMask = np.any(H,axis=-1)
740        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
741        SQfactor = 4.0*SQ*twopisq               #ditto prev.
742        if 'T' in calcControls[hfx+'histType']:
743            if 'P' in calcControls[hfx+'histType']:
744                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
745            else:
746                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
747            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
748            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
749        Uniq = np.inner(H,SGMT)
750        Phi = np.inner(H,SGT)
751        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
752        sinp = np.sin(phase)
753        cosp = np.cos(phase)
754        biso = -SQfactor*Uisodata[:,nxs]
755        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
756        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
757        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
758        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
759        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
760        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
761        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
762            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
763            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
764            if SGData['SGInv']:
765                mphase = np.hstack((phase,-phase))
766            else:
767                mphase = phase
768            mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
769            mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
770            sinm = np.sin(mphase)                               #ditto - match magstrfc.for
771            cosm = np.cos(mphase)                               #ditto
772            HM = np.inner(Bmat.T,H)                             #put into cartesian space
773            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
774            eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
775            Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
776            fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
777            fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
778            fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
779            fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
780            refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
781            refl.T[7] = np.copy(refl.T[9])               
782            refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
783        else:
784            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
785            if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
786                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
787                fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
788            else:
789                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
790                fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
791            fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
792            fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
793            if SGData['SGInv']: #centrosymmetric; B=0
794                fbs[0] *= 0.
795                fas[1] *= 0.
796            if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
797                refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
798                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
799                if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
800                    refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
801            else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
802                if len(TwinLaw) > 1:
803                    refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
804                    refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
805                        np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
806                    refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
807                else:   # checked correct!!
808                    refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
809                    refl.T[7] = np.copy(refl.T[9])               
810                    refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
811#                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
812        iBeg += blkSize
813#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
814   
815def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
816    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
817    faster than StructureFactorDerv - correct for powders/nontwins!!
818    input:
819   
820    :param dict refDict: where
821        'RefList' list where each ref = h,k,l,it,d,...
822        'FF' dict of form factors - filled in below
823    :param np.array G:      reciprocal metric tensor
824    :param str hfx:    histogram id string
825    :param str pfx:    phase id string
826    :param dict SGData: space group info. dictionary output from SpcGroup
827    :param dict calcControls:
828    :param dict parmDict:
829   
830    :returns: dict dFdvDict: dictionary of derivatives
831    '''
832    phfx = pfx.split(':')[0]+hfx
833    ast = np.sqrt(np.diag(G))
834    Mast = twopisq*np.multiply.outer(ast,ast)
835    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
836    SGT = np.array([ops[1] for ops in SGData['SGOps']])
837    FFtables = calcControls['FFtables']
838    BLtables = calcControls['BLtables']
839    Amat,Bmat = G2lat.Gmat2AB(G)
840    nRef = len(refDict['RefList'])
841    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
842        GetAtomFXU(pfx,calcControls,parmDict)
843    if not Xdata.size:          #no atoms in phase!
844        return {}
845    mSize = len(Mdata)
846    FF = np.zeros(len(Tdata))
847    if 'NC' in calcControls[hfx+'histType']:
848        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
849    elif 'X' in calcControls[hfx+'histType']:
850        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
851        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
852    Uij = np.array(G2lat.U6toUij(Uijdata))
853    bij = Mast*Uij.T
854    dFdvDict = {}
855    dFdfr = np.zeros((nRef,mSize))
856    dFdx = np.zeros((nRef,mSize,3))
857    dFdui = np.zeros((nRef,mSize))
858    dFdua = np.zeros((nRef,mSize,6))
859    dFdbab = np.zeros((nRef,2))
860    dFdfl = np.zeros((nRef))
861    Flack = 1.0
862    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
863        Flack = 1.-2.*parmDict[phfx+'Flack']
864#    time0 = time.time()
865#reflection processing begins here - big arrays!
866    iBeg = 0
867    blkSize = 32       #no. of reflections in a block - optimized for speed
868    while iBeg < nRef:
869        iFin = min(iBeg+blkSize,nRef)
870        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
871        H = refl.T[:3].T
872        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
873        SQfactor = 8.0*SQ*np.pi**2
874        if 'T' in calcControls[hfx+'histType']:
875            if 'P' in calcControls[hfx+'histType']:
876                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
877            else:
878                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
879            FP = np.repeat(FP.T,len(SGT),axis=0)
880            FPP = np.repeat(FPP.T,len(SGT),axis=0)
881        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
882        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
883        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
884        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
885        Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
886        Phi = np.inner(H,SGT)
887        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
888        sinp = np.sin(phase)        #refBlk x nOps x nAtoms
889        cosp = np.cos(phase)
890        occ = Mdata*Fdata/len(SGT)
891        biso = -SQfactor*Uisodata[:,nxs]
892        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
893        HbH = np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
894        Tuij = np.where(HbH<1.,np.exp(-HbH),1.0).T
895        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
896        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])      #Nref*Nops,3,3
897        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))     #Nref,Nops,6
898        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
899        if len(FPP.shape) > 1:
900            fotp = np.reshape(FPP,cosp.shape)*Tcorr
901        else:
902            fotp = FPP*Tcorr     
903#            GSASIIpath.IPyBreak()
904        if 'T' in calcControls[hfx+'histType']:
905            fa = np.array([fot*cosp,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
906            fb = np.array([fot*sinp,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
907        else:
908            fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr])
909            fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr])
910        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,refBlk,nTwins)
911        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
912        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,refBlk,nEqv,nAtoms)
913        fbx = np.array([fot*cosp,-fotp*sinp])
914        #sum below is over Uniq
915        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,nAtom) Fdata != 0 avoids /0. problem
916        dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
917        dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
918        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nAtoms)
919        dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
920        # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
921        if not SGData['SGInv']:
922            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
923            dfbdba = np.sum(-sinp*Tcorr,axis=-2)
924            dfadfl = np.sum(np.sum(-fotp*sinp,axis=-1),axis=-1)
925            dfbdfl = np.sum(np.sum(fotp*cosp,axis=-1),axis=-1)
926            dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
927            dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
928            dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2)
929        else:
930            dfbdfr = np.zeros_like(dfadfr)
931            dfbdx = np.zeros_like(dfadx)
932            dfbdui = np.zeros_like(dfadui)
933            dfbdua = np.zeros_like(dfadua)
934            dfbdba = np.zeros_like(dfadba)
935            dfadfl = 0.0
936            dfbdfl = 0.0
937        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
938        SA = fas[0]+fas[1]
939        SB = fbs[0]+fbs[1]
940        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
941            dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/len(SGMT)
942            dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0)
943            dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0)
944            dFdua[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadua+fbs[:,:,nxs,nxs]*dfbdua,axis=0)
945        else:
946            dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
947            dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
948            dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
949            dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])
950            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
951        dFdbab[iBeg:iFin] = 2.*(fas[0,nxs]*np.array([np.sum(dfadba.T*dBabdA,axis=0),np.sum(-dfadba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])+ \
952                            fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
953#        GSASIIpath.IPyBreak()
954        iBeg += blkSize
955#    print 'derv time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
956        #loop over atoms - each dict entry is list of derivatives for all the reflections
957    for i in range(len(Mdata)):
958        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
959        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
960        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
961        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
962        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
963        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
964        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
965        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
966        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
967        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
968        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
969    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
970    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
971    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
972    return dFdvDict
973   
974def StructureFactorDervMag(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
975    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
976    input:
977   
978    :param dict refDict: where
979        'RefList' list where each ref = h,k,l,it,d,...
980        'FF' dict of form factors - filled in below
981    :param np.array G:      reciprocal metric tensor
982    :param str hfx:    histogram id string
983    :param str pfx:    phase id string
984    :param dict SGData: space group info. dictionary output from SpcGroup
985    :param dict calcControls:
986    :param dict parmDict:
987   
988    :returns: dict dFdvDict: dictionary of derivatives
989    '''
990    ast = np.sqrt(np.diag(G))
991    Mast = twopisq*np.multiply.outer(ast,ast)
992    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
993    SGT = np.array([ops[1] for ops in SGData['SGOps']])
994    Ncen = len(SGData['SGCen'])
995    Nops = len(SGMT)*Ncen*(1+SGData['SGInv'])
996    Amat,Bmat = G2lat.Gmat2AB(G)
997    nRef = len(refDict['RefList'])
998    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
999        GetAtomFXU(pfx,calcControls,parmDict)
1000    if not Xdata.size:          #no atoms in phase!
1001        return {}
1002    mSize = len(Mdata)
1003    Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
1004    Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
1005    dGdM = np.repeat(Gdata[:,nxs,:],Nops,axis=1)
1006    Gdata = np.inner(Bmat,Gdata.T)              #convert to crystal space
1007    Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
1008    if SGData['SGInv']:
1009        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
1010    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
1011    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip
1012    Gdata = np.inner(Amat,Gdata.T)              #convert back to cart. space MXYZ, Natoms, NOps
1013    Gdata = np.swapaxes(Gdata,1,2)              # put Natoms last - Mxyz,Nops,Natms
1014    Mag = np.tile(Mag[:,nxs],Nops).#make Mag same length as Gdata
1015    dGdm = (1.-Gdata**2)                        #1/Mag removed - canceled out in dqmx=sum(dqdm*dGdm)
1016    dFdMx = np.zeros((nRef,mSize,3))
1017    Uij = np.array(G2lat.U6toUij(Uijdata))
1018    bij = Mast*Uij.T
1019    dFdvDict = {}
1020    dFdfr = np.zeros((nRef,mSize))
1021    dFdx = np.zeros((nRef,mSize,3))
1022    dFdMx = np.zeros((3,nRef,mSize))
1023    dFdui = np.zeros((nRef,mSize))
1024    dFdua = np.zeros((nRef,mSize,6))
1025    time0 = time.time()
1026#reflection processing begins here - big arrays!
1027    iBeg = 0
1028    blkSize = 32       #no. of reflections in a block - optimized for speed
1029    while iBeg < nRef:
1030        iFin = min(iBeg+blkSize,nRef)
1031        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1032        H = refl.T[:3].T
1033        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1034        SQfactor = 8.0*SQ*np.pi**2
1035        Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
1036        Phi = np.inner(H,SGT)
1037        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1038        occ = Mdata*Fdata/Nops
1039        biso = -SQfactor*Uisodata[:,nxs]
1040        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
1041        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1042        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
1043        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1044        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))
1045        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1046        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1047        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
1048        if SGData['SGInv']:
1049            mphase = np.hstack((phase,-phase))
1050            Uniq = np.hstack((Uniq,-Uniq))      #Nref,Nops,hkl
1051            Hij = np.hstack((Hij,Hij))
1052        else:
1053            mphase = phase
1054        Hij = np.concatenate(np.array([Hij for cen in SGData['SGCen']]),axis=1)
1055        Uniq = np.hstack([Uniq for cen in SGData['SGCen']])
1056        mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
1057        mphase = np.concatenate(mphase,axis=1)              #Nref,Nop,Natm
1058        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
1059        cosm = np.cos(mphase)                               #ditto
1060        HM = np.inner(Bmat.T,H)                             #put into cartesian space
1061        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #unit vector for H
1062        eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
1063        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK
1064        NQ = np.where(np.abs(Q)>0.,1./np.abs(Q),0.)     #this sort of works esp for 1 axis moments
1065#        NQ2 = np.where(np.abs(Q)>0.,1./np.sqrt(np.sum(Q**2,axis=0)),0.)
1066        dqdm = np.array([np.outer(hm,hm)-np.eye(3) for hm in HM.T]).T   #Mxyz,Mxyz,Nref (3x3 matrix)
1067        dqmx = dqdm[:,:,:,nxs,nxs]*dGdm[:,nxs,nxs,:,:]
1068        dqmx2 = np.sum(dqmx,axis=1)   #matrix * vector = vector
1069#        dqmx1 = np.swapaxes(np.swapaxes(np.inner(dqdm.T,dGdm.T),0,1),2,3)
1070        dmx = NQ*Q*dGdM[:,nxs,:,:]-Q*dqmx2                                #*Mag canceled out of dqmx term
1071       
1072        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #Mxyz,Nref,Nop,Natm
1073        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]
1074        fams = np.sum(np.sum(fam,axis=-1),axis=-1)                      #Mxyz,Nref
1075        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)
1076        famx = -Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*sinm[nxs,:,:,:]   #Mxyz,Nref,Nops,Natom
1077        fbmx = Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*cosm[nxs,:,:,:]
1078        #sums below are over Nops - real part
1079        dfadfr = np.sum(fam/occ,axis=2)        #array(Mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem deriv OK
1080        dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*famx[:,:,:,:,nxs],axis=2)          #deriv OK
1081        dfadmx = np.sum(dmx*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:],axis=2)
1082        dfadui = np.sum(-SQfactor[:,nxs,nxs]*fam,axis=2) #array(Ops,refBlk,nAtoms)  deriv OK
1083        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)    #deriv OK? not U12 & U23 in sarc
1084        # imaginary part; array(3,refBlk,nAtom,3) & array(3,refBlk,nAtom,6)
1085        dfbdfr = np.sum(fbm/occ,axis=2)        #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem
1086        dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*fbmx[:,:,:,:,nxs],axis=2)
1087        dfbdmx = np.sum(dmx*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:],axis=2)
1088        dfbdui = np.sum(-SQfactor[:,nxs,nxs]*fbm,axis=2) #array(Ops,refBlk,nAtoms)
1089        dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbm[:,:,:,:,nxs],axis=2)
1090        #accumulate derivatives   
1091        dFdfr[iBeg:iFin] = 2.*np.sum((fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/Nops,axis=0) #ok
1092        dFdx[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadx+fbms[:,:,nxs,nxs]*dfbdx,axis=0)         #ok
1093        dFdMx[:,iBeg:iFin,:] = 2.*(fams[:,:,nxs]*dfadmx+fbms[:,:,nxs]*dfbdmx)                       #problems
1094        dFdui[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui,axis=0)              #ok
1095        dFdua[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua,axis=0)      #problems U12 & U23 in sarc
1096#        GSASIIpath.IPyBreak()
1097        iBeg += blkSize
1098    print ' %d derivative time %.4f\r'%(nRef,time.time()-time0)
1099        #loop over atoms - each dict entry is list of derivatives for all the reflections
1100    for i in range(len(Mdata)):
1101        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1102        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1103        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1104        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1105        dFdvDict[pfx+'AMx:'+str(i)] = dFdMx[0,:,i]
1106        dFdvDict[pfx+'AMy:'+str(i)] = dFdMx[1,:,i]
1107        dFdvDict[pfx+'AMz:'+str(i)] = dFdMx[2,:,i]
1108        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1109        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1110        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1111        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1112        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
1113        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
1114        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
1115#    GSASIIpath.IPyBreak()
1116    return dFdvDict
1117       
1118def StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1119    '''Compute structure factor derivatives on blocks of reflections - for twins only
1120    faster than StructureFactorDervTw
1121    input:
1122   
1123    :param dict refDict: where
1124        'RefList' list where each ref = h,k,l,it,d,...
1125        'FF' dict of form factors - filled in below
1126    :param np.array G:      reciprocal metric tensor
1127    :param str hfx:    histogram id string
1128    :param str pfx:    phase id string
1129    :param dict SGData: space group info. dictionary output from SpcGroup
1130    :param dict calcControls:
1131    :param dict parmDict:
1132   
1133    :returns: dict dFdvDict: dictionary of derivatives
1134    '''
1135    phfx = pfx.split(':')[0]+hfx
1136    ast = np.sqrt(np.diag(G))
1137    Mast = twopisq*np.multiply.outer(ast,ast)
1138    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1139    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1140    FFtables = calcControls['FFtables']
1141    BLtables = calcControls['BLtables']
1142    TwDict = refDict.get('TwDict',{})           
1143    NTL = calcControls[phfx+'NTL']
1144    NM = calcControls[phfx+'TwinNMN']+1
1145    TwinLaw = calcControls[phfx+'TwinLaw']
1146    TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1147    TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1148    nTwin = len(TwinLaw)       
1149    nRef = len(refDict['RefList'])
1150    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1151        GetAtomFXU(pfx,calcControls,parmDict)
1152    if not Xdata.size:          #no atoms in phase!
1153        return {}
1154    mSize = len(Mdata)
1155    FF = np.zeros(len(Tdata))
1156    if 'NC' in calcControls[hfx+'histType']:
1157        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1158    elif 'X' in calcControls[hfx+'histType']:
1159        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1160        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1161    Uij = np.array(G2lat.U6toUij(Uijdata))
1162    bij = Mast*Uij.T
1163    dFdvDict = {}
1164    dFdfr = np.zeros((nRef,nTwin,mSize))
1165    dFdx = np.zeros((nRef,nTwin,mSize,3))
1166    dFdui = np.zeros((nRef,nTwin,mSize))
1167    dFdua = np.zeros((nRef,nTwin,mSize,6))
1168    dFdbab = np.zeros((nRef,nTwin,2))
1169    dFdtw = np.zeros((nRef,nTwin))
1170    time0 = time.time()
1171#reflection processing begins here - big arrays!
1172    iBeg = 0
1173    blkSize = 16       #no. of reflections in a block - optimized for speed
1174    while iBeg < nRef:
1175        iFin = min(iBeg+blkSize,nRef)
1176        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1177        H = refl.T[:3]
1178        H = np.inner(H.T,TwinLaw)   #array(3,nTwins)
1179        TwMask = np.any(H,axis=-1)
1180        for ir in range(blkSize):
1181            iref = ir+iBeg
1182            if iref in TwDict:
1183                for i in TwDict[iref]:
1184                    for n in range(NTL):
1185                        H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1186        TwMask = np.any(H,axis=-1)
1187        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1188        SQfactor = 8.0*SQ*np.pi**2
1189        if 'T' in calcControls[hfx+'histType']:
1190            if 'P' in calcControls[hfx+'histType']:
1191                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1192            else:
1193                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1194            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
1195            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
1196        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1197        Bab = np.repeat(parmDict[phfx+'BabA']*dBabdA,len(SGT)*nTwin)
1198        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1199        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
1200        Uniq = np.inner(H,SGMT)             # (nTwin,nSGOp,3)
1201        Phi = np.inner(H,SGT)
1202        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1203        sinp = np.sin(phase)
1204        cosp = np.cos(phase)
1205        occ = Mdata*Fdata/len(SGT)
1206        biso = -SQfactor*Uisodata[:,nxs]
1207        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
1208        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1209        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1210        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,nTwin,len(SGT),6))
1211        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1212        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*Mdata*Fdata/len(SGMT)
1213        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
1214        fotp = FPP*Tcorr       
1215        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
1216            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(FPP,sinp.shape)*sinp*Tcorr])
1217            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(FPP,cosp.shape)*cosp*Tcorr])
1218        else:
1219            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-FPP*sinp*Tcorr])
1220            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,FPP*cosp*Tcorr])
1221        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
1222        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
1223        if SGData['SGInv']: #centrosymmetric; B=0
1224            fbs[0] *= 0.
1225            fas[1] *= 0.
1226        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,nRef,ntwi,nEqv,nAtoms)
1227        fbx = np.array([fot*cosp,-fotp*sinp])
1228        #sum below is over Uniq
1229        dfadfr = np.sum(np.sum(fa/occ,axis=-2),axis=0)        #array(2,nRef,ntwin,nAtom) Fdata != 0 avoids /0. problem
1230        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1231        dfadui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),axis=0)           
1232        dfadx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fax[:,:,:,:,:,nxs],axis=-3),axis=0) # nRef x nTwin x nAtoms x xyz; sum on ops & A,A'
1233        dfadua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fa[:,:,:,:,:,nxs],axis=-3),axis=0) 
1234        if not SGData['SGInv']:
1235            dfbdfr = np.sum(np.sum(fb/occ,axis=-2),axis=0)        #Fdata != 0 avoids /0. problem
1236            dfadba /= 2.
1237#            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
1238            dfbdui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),axis=0)
1239            dfbdx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fbx[:,:,:,:,:,nxs],axis=-3),axis=0)
1240            dfbdua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fb[:,:,:,:,:,nxs],axis=-3),axis=0)
1241        else:
1242            dfbdfr = np.zeros_like(dfadfr)
1243            dfbdx = np.zeros_like(dfadx)
1244            dfbdui = np.zeros_like(dfadui)
1245            dfbdua = np.zeros_like(dfadua)
1246#            dfbdba = np.zeros_like(dfadba)
1247        SA = fas[0]+fas[1]
1248        SB = fbs[0]+fbs[1]
1249#        GSASIIpath.IPyBreak()
1250        dFdfr[iBeg:iFin] = ((2.*TwMask*SA)[:,:,nxs]*dfadfr+(2.*TwMask*SB)[:,:,nxs]*dfbdfr)*Mdata[nxs,nxs,:]/len(SGMT)
1251        dFdx[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadx+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdx
1252        dFdui[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs]*dfadui+(2.*TwMask*SB)[:,:,nxs]*dfbdui
1253        dFdua[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadua+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdua
1254        if SGData['SGInv']: #centrosymmetric; B=0
1255            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
1256        else:               
1257            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
1258#        dFdbab[iBeg:iFin] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1259#            fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1260        iBeg += blkSize
1261#        GSASIIpath.IPyBreak()
1262    print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)
1263    #loop over atoms - each dict entry is list of derivatives for all the reflections
1264    for i in range(len(Mdata)):     #these all OK
1265        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
1266        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
1267        dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
1268        dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
1269        dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
1270        dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
1271        dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
1272        dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
1273        dFdvDict[pfx+'AU12:'+str(i)] = np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
1274        dFdvDict[pfx+'AU13:'+str(i)] = np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
1275        dFdvDict[pfx+'AU23:'+str(i)] = np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
1276    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1277    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1278    for i in range(nTwin):
1279        dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
1280    return dFdvDict
1281   
1282def SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1283    '''
1284    Compute super structure factors for all h,k,l,m for phase - no twins
1285    puts the result, F^2, in each ref[9] in refList
1286    works on blocks of 32 reflections for speed
1287    input:
1288   
1289    :param dict refDict: where
1290        'RefList' list where each ref = h,k,l,m,it,d,...
1291        'FF' dict of form factors - filed in below
1292    :param np.array G:      reciprocal metric tensor
1293    :param str pfx:    phase id string
1294    :param dict SGData: space group info. dictionary output from SpcGroup
1295    :param dict calcControls:
1296    :param dict ParmDict:
1297
1298    '''
1299    phfx = pfx.split(':')[0]+hfx
1300    ast = np.sqrt(np.diag(G))
1301    Mast = twopisq*np.multiply.outer(ast,ast)   
1302    SGInv = SGData['SGInv']
1303    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1304    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1305    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1306    FFtables = calcControls['FFtables']
1307    BLtables = calcControls['BLtables']
1308    MFtables = calcControls['MFtables']
1309    Flack = 1.0
1310    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1311        Flack = 1.-2.*parmDict[phfx+'Flack']
1312    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1313        GetAtomFXU(pfx,calcControls,parmDict)
1314    if not Xdata.size:          #no atoms in phase!
1315        return
1316    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1317    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1318    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1319    FF = np.zeros(len(Tdata))
1320    if 'NC' in calcControls[hfx+'histType']:
1321        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1322    elif 'X' in calcControls[hfx+'histType']:
1323        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1324        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1325    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1326    bij = Mast*Uij
1327    blkSize = 32       #no. of reflections in a block
1328    nRef = refDict['RefList'].shape[0]
1329    SQ = 1./(2.*refDict['RefList'].T[5])**2
1330    if 'N' in calcControls[hfx+'histType']:
1331        dat = G2el.getBLvalues(BLtables)
1332        refDict['FF']['El'] = dat.keys()
1333        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
1334        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1335        for iel,El in enumerate(refDict['FF']['El']):
1336            if El in MFtables:
1337                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1338    else:
1339        dat = G2el.getFFvalues(FFtables,0.)
1340        refDict['FF']['El'] = dat.keys()
1341        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1342        for iel,El in enumerate(refDict['FF']['El']):
1343            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1344    time0 = time.time()
1345#reflection processing begins here - big arrays!
1346    iBeg = 0
1347    while iBeg < nRef:
1348        iFin = min(iBeg+blkSize,nRef)
1349        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1350        H = refl.T[:4]                          #array(blkSize,4)
1351        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1352        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1353        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1354        Uniq = np.inner(H.T,SSGMT)
1355        UniqP = np.inner(HP.T,SGMT)
1356        Phi = np.inner(H.T,SSGT)
1357        if SGInv:   #if centro - expand HKL sets
1358            Uniq = np.hstack((Uniq,-Uniq))
1359            Phi = np.hstack((Phi,-Phi))
1360            UniqP = np.hstack((UniqP,-UniqP))
1361        if 'T' in calcControls[hfx+'histType']:
1362            if 'P' in calcControls[hfx+'histType']:
1363                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1364            else:
1365                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1366            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
1367            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
1368        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
1369        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1370        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
1371        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1372        sinp = np.sin(phase)
1373        cosp = np.cos(phase)
1374        biso = -SQfactor*Uisodata[:,nxs]
1375        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T
1376        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1377        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1378        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1379        if 'T' in calcControls[hfx+'histType']:
1380            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1381            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1382        else:
1383            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1384            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1385        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1386        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1387        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1388        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1389        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1390#        GSASIIpath.IPyBreak()
1391        if 'P' in calcControls[hfx+'histType']:
1392            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1393#            refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
1394            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1395        else:
1396            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1397            refl.T[8] = np.copy(refl.T[10])               
1398            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1399        iBeg += blkSize
1400    print 'nRef %d time %.4f\r'%(nRef,time.time()-time0)
1401
1402def SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1403    '''
1404    Compute super structure factors for all h,k,l,m for phase - twins only
1405    puts the result, F^2, in each ref[8+im] in refList
1406    works on blocks of 32 reflections for speed
1407    input:
1408   
1409    :param dict refDict: where
1410        'RefList' list where each ref = h,k,l,m,it,d,...
1411        'FF' dict of form factors - filed in below
1412    :param np.array G:      reciprocal metric tensor
1413    :param str pfx:    phase id string
1414    :param dict SGData: space group info. dictionary output from SpcGroup
1415    :param dict calcControls:
1416    :param dict ParmDict:
1417
1418    '''
1419    phfx = pfx.split(':')[0]+hfx
1420    ast = np.sqrt(np.diag(G))
1421    Mast = twopisq*np.multiply.outer(ast,ast)   
1422    SGInv = SGData['SGInv']
1423    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1424    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1425    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1426    FFtables = calcControls['FFtables']
1427    BLtables = calcControls['BLtables']
1428    MFtables = calcControls['MFtables']
1429    Flack = 1.0
1430    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1431        Flack = 1.-2.*parmDict[phfx+'Flack']
1432    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])    #4D?
1433    TwDict = refDict.get('TwDict',{})           
1434    if 'S' in calcControls[hfx+'histType']:
1435        NTL = calcControls[phfx+'NTL']
1436        NM = calcControls[phfx+'TwinNMN']+1
1437        TwinLaw = calcControls[phfx+'TwinLaw']  #this'll have to be 4D also...
1438        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1439        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1440    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1441        GetAtomFXU(pfx,calcControls,parmDict)
1442    if not Xdata.size:          #no atoms in phase!
1443        return
1444    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1445    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1446    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1447    FF = np.zeros(len(Tdata))
1448    if 'NC' in calcControls[hfx+'histType']:
1449        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1450    elif 'X' in calcControls[hfx+'histType']:
1451        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1452        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1453    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1454    bij = Mast*Uij
1455    blkSize = 32       #no. of reflections in a block
1456    nRef = refDict['RefList'].shape[0]
1457    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
1458        SQ = 1./(2.*refDict['RefList'].T[5])**2
1459        if 'N' in calcControls[hfx+'histType']:
1460            dat = G2el.getBLvalues(BLtables)
1461            refDict['FF']['El'] = dat.keys()
1462            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
1463            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1464            for iel,El in enumerate(refDict['FF']['El']):
1465                if El in MFtables:
1466                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1467        else:
1468            dat = G2el.getFFvalues(FFtables,0.)
1469            refDict['FF']['El'] = dat.keys()
1470            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1471            for iel,El in enumerate(refDict['FF']['El']):
1472                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1473    time0 = time.time()
1474#reflection processing begins here - big arrays!
1475    iBeg = 0
1476    while iBeg < nRef:
1477        iFin = min(iBeg+blkSize,nRef)
1478        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1479        H = refl[:,:4]                          #array(blkSize,4)
1480        H3 = refl[:,:3]
1481        HP = H[:,:3]+modQ[nxs,:]*H[:,3:]        #projected hklm to hkl
1482        HP = np.inner(HP,TwinLaw)             #array(blkSize,nTwins,4)
1483        H3 = np.inner(H3,TwinLaw)       
1484        TwMask = np.any(HP,axis=-1)
1485        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1486            for ir in range(blkSize):
1487                iref = ir+iBeg
1488                if iref in TwDict:
1489                    for i in TwDict[iref]:
1490                        for n in range(NTL):
1491                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1492                            H3[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1493            TwMask = np.any(HP,axis=-1)
1494        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1495        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1496        Uniq = np.inner(H,SSGMT)
1497        Uniq3 = np.inner(H3,SGMT)
1498        UniqP = np.inner(HP,SGMT)
1499        Phi = np.inner(H,SSGT)
1500        if SGInv:   #if centro - expand HKL sets
1501            Uniq = np.hstack((Uniq,-Uniq))
1502            Uniq3 = np.hstack((Uniq3,-Uniq3))
1503            Phi = np.hstack((Phi,-Phi))
1504            UniqP = np.hstack((UniqP,-UniqP))
1505        if 'T' in calcControls[hfx+'histType']:
1506            if 'P' in calcControls[hfx+'histType']:
1507                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1508            else:
1509                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1510            FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1511            FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1512        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw))
1513        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1514        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1515        phase = twopi*(np.inner(Uniq3,(dXdata.T+Xdata.T))-Phi[:,nxs,:,nxs])
1516        sinp = np.sin(phase)
1517        cosp = np.cos(phase)
1518        biso = -SQfactor*Uisodata[:,nxs]
1519        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1520        HbH = -np.sum(UniqP[:,:,:,nxs]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1521        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1522        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1523#        GSASIIpath.IPyBreak()
1524        if 'T' in calcControls[hfx+'histType']:
1525            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1526            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1527        else:
1528            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1529            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1530        GfpuA = G2mth.ModulationTw(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1531        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1532        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1533        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1534        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1535        refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
1536        refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
1537            np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
1538        refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
1539        iBeg += blkSize
1540    print 'nRef %d time %.4f\r'%(nRef,time.time()-time0)
1541
1542def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1543    '''
1544    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
1545    input:
1546   
1547    :param dict refDict: where
1548        'RefList' list where each ref = h,k,l,m,it,d,...
1549        'FF' dict of form factors - filled in below
1550    :param int im: = 1 (could be eliminated)
1551    :param np.array G:      reciprocal metric tensor
1552    :param str hfx:    histogram id string
1553    :param str pfx:    phase id string
1554    :param dict SGData: space group info. dictionary output from SpcGroup
1555    :param dict SSGData: super space group info.
1556    :param dict calcControls:
1557    :param dict ParmDict:
1558   
1559    :returns: dict dFdvDict: dictionary of derivatives
1560    '''
1561    phfx = pfx.split(':')[0]+hfx
1562    ast = np.sqrt(np.diag(G))
1563    Mast = twopisq*np.multiply.outer(ast,ast)
1564    SGInv = SGData['SGInv']
1565    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1566    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1567    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1568    FFtables = calcControls['FFtables']
1569    BLtables = calcControls['BLtables']
1570    nRef = len(refDict['RefList'])
1571    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1572        GetAtomFXU(pfx,calcControls,parmDict)
1573    if not Xdata.size:          #no atoms in phase!
1574        return {}
1575    mSize = len(Mdata)  #no. atoms
1576    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1577    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1578    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
1579    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1580    FF = np.zeros(len(Tdata))
1581    if 'NC' in calcControls[hfx+'histType']:
1582        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1583    elif 'X' in calcControls[hfx+'histType']:
1584        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1585        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1586    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1587    bij = Mast*Uij
1588    if not len(refDict['FF']):
1589        if 'N' in calcControls[hfx+'histType']:
1590            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1591        else:
1592            dat = G2el.getFFvalues(FFtables,0.)       
1593        refDict['FF']['El'] = dat.keys()
1594        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1595    dFdvDict = {}
1596    dFdfr = np.zeros((nRef,mSize))
1597    dFdx = np.zeros((nRef,mSize,3))
1598    dFdui = np.zeros((nRef,mSize))
1599    dFdua = np.zeros((nRef,mSize,6))
1600    dFdbab = np.zeros((nRef,2))
1601    dFdfl = np.zeros((nRef))
1602    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1603    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1604    dFdGz = np.zeros((nRef,mSize,5))
1605    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1606    Flack = 1.0
1607    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1608        Flack = 1.-2.*parmDict[phfx+'Flack']
1609    time0 = time.time()
1610    nRef = len(refDict['RefList'])/100
1611    for iref,refl in enumerate(refDict['RefList']):
1612        if 'T' in calcControls[hfx+'histType']:
1613            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1614        H = np.array(refl[:4])
1615        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
1616        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1617        SQfactor = 8.0*SQ*np.pi**2
1618        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1619        Bab = parmDict[phfx+'BabA']*dBabdA
1620        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1621        FF = refDict['FF']['FF'][iref].T[Tindx]
1622        Uniq = np.inner(H,SSGMT)
1623        Phi = np.inner(H,SSGT)
1624        UniqP = np.inner(HP,SGMT)
1625        if SGInv:   #if centro - expand HKL sets
1626            Uniq = np.vstack((Uniq,-Uniq))
1627            Phi = np.hstack((Phi,-Phi))
1628            UniqP = np.vstack((UniqP,-UniqP))
1629        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1630        sinp = np.sin(phase)
1631        cosp = np.cos(phase)
1632        occ = Mdata*Fdata/Uniq.shape[0]
1633        biso = -SQfactor*Uisodata[:,nxs]
1634        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0],axis=1).T    #ops x atoms
1635        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
1636        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1637        Hij = np.array([G2lat.UijtoU6(uij) for uij in Hij])                     #atoms x 6
1638        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1639        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1640        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
1641        fotp = FPP*Tcorr            #ops x atoms
1642        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1643        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1644        # GfpuA is 2 x ops x atoms
1645        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1646        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
1647        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1648        fag = fa*GfpuA[0]-fb*GfpuA[1]
1649        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1650       
1651        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
1652        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1653        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
1654        fbx = np.array([fot*cosp,-fotp*sinp])
1655        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1656        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1657        #sum below is over Uniq
1658        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
1659        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
1660        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1661        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1662        dfadui = np.sum(-SQfactor*fag,axis=1)
1663        dfbdui = np.sum(-SQfactor*fbg,axis=1)
1664        dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
1665        dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
1666        dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
1667        dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
1668        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1669        dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1670        dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1671        dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1672        dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1673        dfadGz = np.sum(fa[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]-fb[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
1674        dfbdGz = np.sum(fb[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]+fa[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
1675        dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
1676        dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
1677        if not SGData['SGInv']:   #Flack derivative
1678            dfadfl = np.sum(-FPP*Tcorr*sinp)
1679            dfbdfl = np.sum(FPP*Tcorr*cosp)
1680        else:
1681            dfadfl = 1.0
1682            dfbdfl = 1.0
1683#        GSASIIpath.IPyBreak()
1684        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
1685        SA = fas[0]+fas[1]      #float = A+A'
1686        SB = fbs[0]+fbs[1]      #float = B+B'
1687        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
1688            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1689            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
1690                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1691            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
1692                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1693            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
1694                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1695            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
1696                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1697            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
1698                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
1699            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
1700                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
1701            dFdGz[iref] = 2.*(fas[0]*dfadGz[0]+fas[1]*dfadGz[1])+  \
1702                2.*(fbs[0]*dfbdGz[0]+fbs[1]*dfbdGz[1])
1703            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
1704                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
1705        else:                       #OK, I think
1706            dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
1707            dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
1708            dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
1709            dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
1710            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1711                           
1712            dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
1713            dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
1714            dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1])      #array(nRef,natom,5)
1715            dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
1716#            GSASIIpath.IPyBreak()
1717        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1718            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1719        #loop over atoms - each dict entry is list of derivatives for all the reflections
1720        if not iref%100 :
1721            print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
1722    for i in range(len(Mdata)):     #loop over atoms
1723        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1724        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1725        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1726        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1727        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1728        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1729        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1730        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1731        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
1732        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
1733        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
1734        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
1735            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
1736            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
1737        nx = 0
1738        if waveTypes[i] in ['Block','ZigZag']:
1739            nx = 1 
1740            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
1741            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
1742            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
1743            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
1744            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
1745        for j in range(XSSdata.shape[1]-nx):       #loop over waves
1746            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
1747            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
1748            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
1749            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
1750            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
1751            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
1752        for j in range(USSdata.shape[1]):       #loop over waves
1753            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
1754            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
1755            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
1756            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
1757            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
1758            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
1759            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
1760            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
1761            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
1762            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
1763            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
1764            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
1765           
1766#        GSASIIpath.IPyBreak()
1767    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
1768    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1769    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1770    return dFdvDict
1771
1772def SStructureFactorDerv2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1773    'Needs a doc string - no twins'
1774    phfx = pfx.split(':')[0]+hfx
1775    ast = np.sqrt(np.diag(G))
1776    Mast = twopisq*np.multiply.outer(ast,ast)
1777    SGInv = SGData['SGInv']
1778    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1779    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1780    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1781    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1782    FFtables = calcControls['FFtables']
1783    BLtables = calcControls['BLtables']
1784    nRef = len(refDict['RefList'])
1785    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1786        GetAtomFXU(pfx,calcControls,parmDict)
1787    if not Xdata.size:          #no atoms in phase!
1788        return {}
1789    mSize = len(Mdata)  #no. atoms
1790    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1791    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1792    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
1793    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1794    FF = np.zeros(len(Tdata))
1795    if 'NC' in calcControls[hfx+'histType']:
1796        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1797    elif 'X' in calcControls[hfx+'histType']:
1798        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1799        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1800    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1801    bij = Mast*Uij
1802    if not len(refDict['FF']):
1803        if 'N' in calcControls[hfx+'histType']:
1804            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1805        else:
1806            dat = G2el.getFFvalues(FFtables,0.)       
1807        refDict['FF']['El'] = dat.keys()
1808        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1809    dFdvDict = {}
1810    dFdfr = np.zeros((nRef,mSize))
1811    dFdx = np.zeros((nRef,mSize,3))
1812    dFdui = np.zeros((nRef,mSize))
1813    dFdua = np.zeros((nRef,mSize,6))
1814    dFdbab = np.zeros((nRef,2))
1815    dFdfl = np.zeros((nRef))
1816    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1817    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1818    dFdGz = np.zeros((nRef,mSize,5))
1819    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1820    Flack = 1.0
1821    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1822        Flack = 1.-2.*parmDict[phfx+'Flack']
1823    time0 = time.time()
1824    iBeg = 0
1825    blkSize = 4       #no. of reflections in a block - optimized for speed
1826    while iBeg < nRef:
1827        iFin = min(iBeg+blkSize,nRef)
1828        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1829        H = refl.T[:4]
1830        HP = H[:3].T+modQ*H.T[:,3:]            #projected hklm to hkl
1831        SQ = 1./(2.*refl.T[4+im])**2             # or (sin(theta)/lambda)**2
1832        SQfactor = 8.0*SQ*np.pi**2
1833        if 'T' in calcControls[hfx+'histType']:
1834            if 'P' in calcControls[hfx+'histType']:
1835                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[15])
1836            else:
1837                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[13])
1838            FP = np.repeat(FP.T,len(SGT),axis=0)
1839            FPP = np.repeat(FPP.T,len(SGT),axis=0)
1840#        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1841        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
1842        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1843        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
1844        Uniq = np.inner(H.T,SSGMT)
1845        Phi = np.inner(H.T,SSGT)
1846        UniqP = np.inner(HP,SGMT)
1847        if SGInv:   #if centro - expand HKL sets
1848            Uniq = np.hstack((Uniq,-Uniq))
1849            Phi = np.hstack((Phi,-Phi))
1850            UniqP = np.hstack((UniqP,-UniqP))
1851            FF = np.vstack((FF,FF))
1852            Bab = np.concatenate((Bab,Bab))
1853        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata+Xdata).T)+Phi[:,:,nxs])
1854        sinp = np.sin(phase)
1855        cosp = np.cos(phase)
1856        occ = Mdata*Fdata/Uniq.shape[1]
1857        biso = -SQfactor*Uisodata[:,nxs]
1858        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T    #ops x atoms
1859        HbH = -np.sum(UniqP[:,:,nxs,:3]*np.inner(UniqP[:,:,:3],bij),axis=-1)  #ops x atoms
1860        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in np.reshape(UniqP,(-1,3))]) #atoms x 3x3
1861        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(iFin-iBeg,-1,6))                     #atoms x 6
1862        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1863#        GSASIIpath.IPyBreak()
1864        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1865        fot = np.reshape(FF+FP[nxs,:]-Bab[:,nxs],cosp.shape)*Tcorr     #ops x atoms
1866        fotp = FPP*Tcorr            #ops x atoms
1867        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1868        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv2(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1869        # GfpuA is 2 x ops x atoms
1870        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1871        fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
1872        fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1873        fag = fa*GfpuA[0]-fb*GfpuA[1]
1874        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1875       
1876        fas = np.sum(np.sum(fag,axis=-1),axis=-1)     # 2 x refBlk
1877        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1878        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
1879        fbx = np.array([fot*cosp,-fotp*sinp])
1880        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1881        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1882        #sum below is over Uniq
1883        dfadfr = np.sum(fag/occ,axis=-2)        #Fdata != 0 ever avoids /0. problem
1884        dfbdfr = np.sum(fbg/occ,axis=-2)        #Fdata != 0 avoids /0. problem
1885#        dfadba = np.sum(-cosp*Tcorr,axis=-2)
1886#        dfbdba = np.sum(-sinp*Tcorr,axis=-2)
1887        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fag,axis=-2)
1888        dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fbg,axis=-2)
1889        dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:3]*fax[:,:,:,:,nxs],axis=-3)  #2 refBlk x x nAtom x 3xyz; sum nOps
1890        dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:3]*fbx[:,:,:,:,nxs],axis=-3)  #2 refBlk x x nAtom x 3xyz; sum nOps
1891        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fag[:,:,:,:,nxs],axis=-3)             #2 x nAtom x 6Uij; sum nOps
1892        dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbg[:,:,:,:,nxs],axis=-3)             #2 x nAtom x 6Uij; sum nOps
1893        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1894        dfadGf = np.sum(fa[:,:,:,:,nxs,nxs]*dGdf[0][nxs,:,nxs,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdf[1][nxs,:,nxs,:,:,:],axis=2)
1895        dfbdGf = np.sum(fb[:,:,:,:,nxs,nxs]*dGdf[0][nxs,:,nxs,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdf[1][nxs,:,nxs,:,:,:],axis=2)
1896        dfadGx = np.sum(fa[:,:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:,:],axis=2)
1897        dfbdGx = np.sum(fb[:,:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:,:],axis=2)
1898        dfadGz = np.sum(fa[:,:,:,:,nxs]*dGdz[0][nxs,:,:,:,:]-fb[:,:,:,:,nxs]*dGdz[1][nxs,:,:,:,:],axis=2)
1899        dfbdGz = np.sum(fb[:,:,:,:,nxs]*dGdz[0][nxs,:,:,:,:]+fa[:,:,:,:,nxs]*dGdz[1][nxs,:,:,:,:],axis=2)
1900        dfadGu = np.sum(fa[:,:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=2)
1901        dfbdGu = np.sum(fb[:,:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=2)   
1902        if not SGData['SGInv']:   #Flack derivative
1903            dfadfl = np.sum(np.sum(-FPP*Tcorr*sinp,axis=-1),axis=-1)
1904            dfbdfl = np.sum(np.sum(FPP*Tcorr*cosp,axis=-1),axis=-1)
1905        else:
1906            dfadfl = 1.0
1907            dfbdfl = 1.0
1908        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
1909        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
1910        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
1911        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
1912            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1913            dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(Uniq)+   \
1914                2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(Uniq)
1915            dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+  \
1916                2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1])
1917            dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+   \
1918                2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1])
1919            dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+   \
1920                2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1])
1921            dFdGf[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGf[0]+fas[1,:,nxs,nxs,nxs]*dfadGf[1])+  \
1922                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGf[0]+fbs[1,:,nxs,nxs,nxs]*dfbdGf[1])
1923            dFdGx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGx[0]+fas[1,:,nxs,nxs,nxs]*dfadGx[1])+  \
1924                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGx[0]-fbs[1,:,nxs,nxs,nxs]*dfbdGx[1])
1925            dFdGz[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadGz[0]+fas[1,:,nxs,nxs]*dfadGz[1])+  \
1926                2.*(fbs[0,:,nxs,nxs]*dfbdGz[0]+fbs[1,:,nxs,nxs]*dfbdGz[1])
1927            dFdGu[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGu[0]+fas[1,:,nxs,nxs,nxs]*dfadGu[1])+  \
1928                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGu[0]+fbs[1,:,nxs,nxs,nxs]*dfbdGu[1])
1929        else:                       #OK, I think
1930            dFdfr[iBeg:iFin] = 2.*(SA[:,nxs]*(dfadfr[0]+dfadfr[1])+SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq) #array(nRef,nAtom)
1931            dFdx[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1]))    #array(nRef,nAtom,3)
1932            dFdui[iBeg:iFin] = 2.*(SA[:,nxs]*(dfadui[0]+dfadui[1])+SB[:,nxs]*(dfbdui[0]+dfbdui[1]))   #array(nRef,nAtom)
1933            dFdua[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1]))    #array(nRef,nAtom,6)
1934            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1935                           
1936            dFdGf[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGf[0]+SB[:,nxs,nxs,nxs]*dfbdGf[1])      #array(nRef,natom,nwave,2)
1937            dFdGx[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGx[0]+SB[:,nxs,nxs,nxs]*dfbdGx[1])      #array(nRef,natom,nwave,6)
1938            dFdGz[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*dfadGz[0]+SB[:,nxs,nxs]*dfbdGz[1])      #array(nRef,natom,5)
1939            dFdGu[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGu[0]+SB[:,nxs,nxs,nxs]*dfbdGu[1])      #array(nRef,natom,nwave,12)
1940#            GSASIIpath.IPyBreak()
1941#        dFdbab[iBeg:iFin] = 2.*fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1942#            2.*fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1943        #loop over atoms - each dict entry is list of derivatives for all the reflections
1944        print ' %d derivative time %.4f\r'%(iBeg,time.time()-time0),
1945        iBeg += blkSize
1946    for i in range(len(Mdata)):     #loop over atoms
1947        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1948        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1949        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1950        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1951        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1952        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1953        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1954        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1955        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
1956        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
1957        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
1958        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
1959            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
1960            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
1961        nx = 0
1962        if waveTypes[i] in ['Block','ZigZag']:
1963            nx = 1 
1964            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
1965            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
1966            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
1967            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
1968            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
1969        for j in range(XSSdata.shape[1]-nx):       #loop over waves
1970            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
1971            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
1972            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
1973            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
1974            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
1975            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
1976        for j in range(USSdata.shape[1]):       #loop over waves
1977            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
1978            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
1979            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
1980            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
1981            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
1982            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
1983            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
1984            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
1985            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
1986            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
1987            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
1988            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
1989           
1990#        GSASIIpath.IPyBreak()
1991    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1992    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1993    return dFdvDict
1994   
1995def SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1996    'Needs a doc string'
1997    phfx = pfx.split(':')[0]+hfx
1998    ast = np.sqrt(np.diag(G))
1999    Mast = twopisq*np.multiply.outer(ast,ast)
2000    SGInv = SGData['SGInv']
2001    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2002    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2003    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2004    FFtables = calcControls['FFtables']
2005    BLtables = calcControls['BLtables']
2006    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
2007    TwDict = refDict.get('TwDict',{})           
2008    if 'S' in calcControls[hfx+'histType']:
2009        NTL = calcControls[phfx+'NTL']
2010        NM = calcControls[phfx+'TwinNMN']+1
2011        TwinLaw = calcControls[phfx+'TwinLaw']
2012        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
2013    nTwin = len(TwinLaw)       
2014    nRef = len(refDict['RefList'])
2015    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
2016        GetAtomFXU(pfx,calcControls,parmDict)
2017    if not Xdata.size:          #no atoms in phase!
2018        return {}
2019    mSize = len(Mdata)  #no. atoms
2020    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
2021    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
2022    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
2023    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2024    FF = np.zeros(len(Tdata))
2025    if 'NC' in calcControls[hfx+'histType']:
2026        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
2027    elif 'X' in calcControls[hfx+'histType']:
2028        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2029        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2030    Uij = np.array(G2lat.U6toUij(Uijdata)).T
2031    bij = Mast*Uij
2032    if not len(refDict['FF']):
2033        if 'N' in calcControls[hfx+'histType']:
2034            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
2035        else:
2036            dat = G2el.getFFvalues(FFtables,0.)       
2037        refDict['FF']['El'] = dat.keys()
2038        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
2039    dFdvDict = {}
2040    dFdfr = np.zeros((nRef,nTwin,mSize))
2041    dFdx = np.zeros((nRef,nTwin,mSize,3))
2042    dFdui = np.zeros((nRef,nTwin,mSize))
2043    dFdua = np.zeros((nRef,nTwin,mSize,6))
2044    dFdbab = np.zeros((nRef,nTwin,2))
2045    dFdtw = np.zeros((nRef,nTwin))
2046    dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
2047    dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
2048    dFdGz = np.zeros((nRef,nTwin,mSize,5))
2049    dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
2050    Flack = 1.0
2051    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
2052        Flack = 1.-2.*parmDict[phfx+'Flack']
2053    time0 = time.time()
2054    nRef = len(refDict['RefList'])/100
2055    for iref,refl in enumerate(refDict['RefList']):
2056        if 'T' in calcControls[hfx+'histType']:
2057            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
2058        H = np.array(refl[:4])
2059        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
2060        H = np.inner(H.T,TwinLaw)   #maybe array(4,nTwins) or (4)
2061        TwMask = np.any(H,axis=-1)
2062        if TwinLaw.shape[0] > 1 and TwDict:
2063            if iref in TwDict:
2064                for i in TwDict[iref]:
2065                    for n in range(NTL):
2066                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
2067            TwMask = np.any(H,axis=-1)
2068        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
2069        SQfactor = 8.0*SQ*np.pi**2
2070        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2071        Bab = parmDict[phfx+'BabA']*dBabdA
2072        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
2073        FF = refDict['FF']['FF'][iref].T[Tindx]
2074        Uniq = np.inner(H,SSGMT)
2075        Phi = np.inner(H,SSGT)
2076        UniqP = np.inner(HP,SGMT)
2077        if SGInv:   #if centro - expand HKL sets
2078            Uniq = np.vstack((Uniq,-Uniq))
2079            Phi = np.hstack((Phi,-Phi))
2080            UniqP = np.vstack((UniqP,-UniqP))
2081        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
2082        sinp = np.sin(phase)
2083        cosp = np.cos(phase)
2084        occ = Mdata*Fdata/Uniq.shape[0]
2085        biso = -SQfactor*Uisodata[:,nxs]
2086        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
2087        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
2088        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
2089        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(nTwin,-1,6)))
2090        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
2091        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
2092        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
2093        fotp = FPP*Tcorr            #ops x atoms
2094        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
2095        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
2096        # GfpuA is 2 x ops x atoms
2097        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
2098        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
2099        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
2100        fag = fa*GfpuA[0]-fb*GfpuA[1]
2101        fbg = fb*GfpuA[0]+fa*GfpuA[1]
2102       
2103        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
2104        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
2105        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x twin x ops x atoms
2106        fbx = np.array([fot*cosp,-fotp*sinp])
2107        fax = fax*GfpuA[0]-fbx*GfpuA[1]
2108        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
2109        #sum below is over Uniq
2110        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
2111        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
2112        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
2113        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
2114        dfadui = np.sum(-SQfactor*fag,axis=1)
2115        dfbdui = np.sum(-SQfactor*fbg,axis=1)
2116        dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2117        dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
2118        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2119        dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2120        # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
2121        dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2122        dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2123        dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2124        dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2125        dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2126        dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2127        dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2128        dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2129#        GSASIIpath.IPyBreak()
2130        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
2131        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
2132        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
2133        dFdfr[iref] = [2.*TwMask[it]*(SA[it]*dfadfr[0][it]+SA[it]*dfadfr[1][it]+SB[it]*dfbdfr[0][it]+SB[it]*dfbdfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
2134        dFdx[iref] = [2.*TwMask[it]*(SA[it]*dfadx[it][0]+SA[it]*dfadx[it][1]+SB[it]*dfbdx[it][0]+SB[it]*dfbdx[it][1]) for it in range(nTwin)]
2135        dFdui[iref] = [2.*TwMask[it]*(SA[it]*dfadui[it][0]+SA[it]*dfadui[it][1]+SB[it]*dfbdui[it][0]+SB[it]*dfbdui[it][1]) for it in range(nTwin)]
2136        dFdua[iref] = [2.*TwMask[it]*(SA[it]*dfadua[it][0]+SA[it]*dfadua[it][1]+SB[it]*dfbdua[it][0]+SB[it]*dfbdua[it][1]) for it in range(nTwin)]
2137        dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
2138
2139        dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
2140        dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
2141        dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
2142        dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
2143#            GSASIIpath.IPyBreak()
2144        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
2145            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2146        #loop over atoms - each dict entry is list of derivatives for all the reflections
2147        if not iref%100 :
2148            print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
2149    for i in range(len(Mdata)):     #loop over atoms
2150        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2151        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2152        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2153        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2154        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2155        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2156        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2157        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2158        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
2159        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
2160        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
2161        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2162            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2163            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2164        nx = 0
2165        if waveTypes[i] in ['Block','ZigZag']:
2166            nx = 1 
2167            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
2168            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
2169            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
2170            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
2171            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
2172        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2173            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2174            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2175            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2176            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2177            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2178            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2179        for j in range(USSdata.shape[1]):       #loop over waves
2180            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2181            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2182            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2183            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
2184            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
2185            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
2186            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2187            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2188            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2189            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
2190            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
2191            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
2192           
2193#        GSASIIpath.IPyBreak()
2194    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2195    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2196    return dFdvDict
2197   
2198def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
2199    ''' Single crystal extinction function; returns extinction & derivative
2200    '''
2201    extCor = 1.0
2202    dervDict = {}
2203    dervCor = 1.0
2204    if calcControls[phfx+'EType'] != 'None':
2205        SQ = 1/(4.*ref[4+im]**2)
2206        if 'C' in parmDict[hfx+'Type']:           
2207            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
2208        else:   #'T'
2209            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
2210        if 'SXC' in parmDict[hfx+'Type']:
2211            AV = 7.9406e5/parmDict[pfx+'Vol']**2
2212            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2213            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
2214            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
2215        elif 'SNT' in parmDict[hfx+'Type']:
2216            AV = 1.e7/parmDict[pfx+'Vol']**2
2217            PL = SQ
2218            PLZ = AV*ref[9+im]*ref[12+im]**2
2219        elif 'SNC' in parmDict[hfx+'Type']:
2220            AV = 1.e7/parmDict[pfx+'Vol']**2
2221            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2222            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
2223           
2224        if 'Primary' in calcControls[phfx+'EType']:
2225            PLZ *= 1.5
2226        else:
2227            if 'C' in parmDict[hfx+'Type']:
2228                PLZ *= calcControls[phfx+'Tbar']
2229            else: #'T'
2230                PLZ *= ref[13+im]      #t-bar
2231        if 'Primary' in calcControls[phfx+'EType']:
2232            PLZ *= 1.5
2233            PSIG = parmDict[phfx+'Ep']
2234        elif 'I & II' in calcControls[phfx+'EType']:
2235            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
2236        elif 'Type II' in calcControls[phfx+'EType']:
2237            PSIG = parmDict[phfx+'Es']
2238        else:       # 'Secondary Type I'
2239            PSIG = parmDict[phfx+'Eg']/PL
2240           
2241        AG = 0.58+0.48*cos2T+0.24*cos2T**2
2242        AL = 0.025+0.285*cos2T
2243        BG = 0.02-0.025*cos2T
2244        BL = 0.15-0.2*(0.75-cos2T)**2
2245        if cos2T < 0.:
2246            BL = -0.45*cos2T
2247        CG = 2.
2248        CL = 2.
2249        PF = PLZ*PSIG
2250       
2251        if 'Gaussian' in calcControls[phfx+'EApprox']:
2252            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
2253            extCor = np.sqrt(PF4)
2254            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
2255        else:
2256            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
2257            extCor = np.sqrt(PF4)
2258            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
2259
2260        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
2261        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
2262            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
2263        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
2264            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
2265        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
2266            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
2267               
2268    return 1./extCor,dervDict,dervCor
2269   
2270def Dict2Values(parmdict, varylist):
2271    '''Use before call to leastsq to setup list of values for the parameters
2272    in parmdict, as selected by key in varylist'''
2273    return [parmdict[key] for key in varylist] 
2274   
2275def Values2Dict(parmdict, varylist, values):
2276    ''' Use after call to leastsq to update the parameter dictionary with
2277    values corresponding to keys in varylist'''
2278    parmdict.update(zip(varylist,values))
2279   
2280def GetNewCellParms(parmDict,varyList):
2281    'Needs a doc string'
2282    newCellDict = {}
2283    Anames = ['A'+str(i) for i in range(6)]
2284    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
2285    for item in varyList:
2286        keys = item.split(':')
2287        if keys[2] in Ddict:
2288            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
2289            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
2290            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2291    return newCellDict          # is e.g. {'0::D11':A0-D11}
2292   
2293def ApplyXYZshifts(parmDict,varyList):
2294    '''
2295    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2296   
2297    :param dict parmDict: parameter dictionary
2298    :param list varyList: list of variables (not used!)
2299    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
2300
2301    '''
2302    newAtomDict = {}
2303    for item in parmDict:
2304        if 'dA' in item:
2305            parm = ''.join(item.split('d'))
2306            parmDict[parm] += parmDict[item]
2307            newAtomDict[item] = [parm,parmDict[parm]]
2308    return newAtomDict
2309   
2310def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2311    'Spherical harmonics texture'
2312    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2313    if 'T' in calcControls[hfx+'histType']:
2314        tth = parmDict[hfx+'2-theta']
2315    else:
2316        tth = refl[5+im]
2317    odfCor = 1.0
2318    H = refl[:3]
2319    cell = G2lat.Gmat2cell(g)
2320    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2321    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2322    phi,beta = G2lat.CrsAng(H,cell,SGData)
2323    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2324    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2325    for item in SHnames:
2326        L,M,N = eval(item.strip('C'))
2327        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2328        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2329        Lnorm = G2lat.Lnorm(L)
2330        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2331    return odfCor
2332   
2333def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2334    'Spherical harmonics texture derivatives'
2335    if 'T' in calcControls[hfx+'histType']:
2336        tth = parmDict[hfx+'2-theta']
2337    else:
2338        tth = refl[5+im]
2339    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2340    odfCor = 1.0
2341    dFdODF = {}
2342    dFdSA = [0,0,0]
2343    H = refl[:3]
2344    cell = G2lat.Gmat2cell(g)
2345    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2346    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2347    phi,beta = G2lat.CrsAng(H,cell,SGData)
2348    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
2349    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2350    for item in SHnames:
2351        L,M,N = eval(item.strip('C'))
2352        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2353        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2354        Lnorm = G2lat.Lnorm(L)
2355        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2356        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2357        for i in range(3):
2358            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2359    return odfCor,dFdODF,dFdSA
2360   
2361def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2362    'spherical harmonics preferred orientation (cylindrical symmetry only)'
2363    if 'T' in calcControls[hfx+'histType']:
2364        tth = parmDict[hfx+'2-theta']
2365    else:
2366        tth = refl[5+im]
2367    odfCor = 1.0
2368    H = refl[:3]
2369    cell = G2lat.Gmat2cell(g)
2370    Sangls = [0.,0.,0.]
2371    if 'Bragg' in calcControls[hfx+'instType']:
2372        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2373        IFCoup = True
2374    else:
2375        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2376        IFCoup = False
2377    phi,beta = G2lat.CrsAng(H,cell,SGData)
2378    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2379    SHnames = calcControls[phfx+'SHnames']
2380    for item in SHnames:
2381        L,N = eval(item.strip('C'))
2382        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2383        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2384        Lnorm = G2lat.Lnorm(L)
2385        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2386    return np.squeeze(odfCor)
2387   
2388def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2389    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
2390    if 'T' in calcControls[hfx+'histType']:
2391        tth = parmDict[hfx+'2-theta']
2392    else:
2393        tth = refl[5+im]
2394    odfCor = 1.0
2395    dFdODF = {}
2396    H = refl[:3]
2397    cell = G2lat.Gmat2cell(g)
2398    Sangls = [0.,0.,0.]
2399    if 'Bragg' in calcControls[hfx+'instType']:
2400        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2401        IFCoup = True
2402    else:
2403        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2404        IFCoup = False
2405    phi,beta = G2lat.CrsAng(H,cell,SGData)
2406    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2407    SHnames = calcControls[phfx+'SHnames']
2408    for item in SHnames:
2409        L,N = eval(item.strip('C'))
2410        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2411        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2412        Lnorm = G2lat.Lnorm(L)
2413        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2414        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
2415    return odfCor,dFdODF
2416   
2417def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2418    'March-Dollase preferred orientation correction'
2419    POcorr = 1.0
2420    MD = parmDict[phfx+'MD']
2421    if MD != 1.0:
2422        MDAxis = calcControls[phfx+'MDAxis']
2423        sumMD = 0
2424        for H in uniq:           
2425            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2426            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2427            sumMD += A**3
2428        POcorr = sumMD/len(uniq)
2429    return POcorr
2430   
2431def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2432    'Needs a doc string'
2433    POcorr = 1.0
2434    POderv = {}
2435    if calcControls[phfx+'poType'] == 'MD':
2436        MD = parmDict[phfx+'MD']
2437        MDAxis = calcControls[phfx+'MDAxis']
2438        sumMD = 0
2439        sumdMD = 0
2440        for H in uniq:           
2441            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2442            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2443            sumMD += A**3
2444            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2445        POcorr = sumMD/len(uniq)
2446        POderv[phfx+'MD'] = sumdMD/len(uniq)
2447    else:   #spherical harmonics
2448        if calcControls[phfx+'SHord']:
2449            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2450    return POcorr,POderv
2451   
2452def GetAbsorb(refl,im,hfx,calcControls,parmDict):
2453    'Needs a doc string'
2454    if 'Debye' in calcControls[hfx+'instType']:
2455        if 'T' in calcControls[hfx+'histType']:
2456            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2457        else:
2458            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2459    else:
2460        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
2461   
2462def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
2463    'Needs a doc string'
2464    if 'Debye' in calcControls[hfx+'instType']:
2465        if 'T' in calcControls[hfx+'histType']:
2466            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2467        else:
2468            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2469    else:
2470        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
2471       
2472def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2473    'Needs a doc string'
2474    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2475    pi2 = np.sqrt(2./np.pi)
2476    if 'T' in calcControls[hfx+'histType']:
2477        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2478        wave = refl[14+im]
2479    else:   #'C'W
2480        sth2 = sind(refl[5+im]/2.)**2
2481        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2482    c2th = 1.-2.0*sth2
2483    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2484    if 'X' in calcControls[hfx+'histType']:
2485        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2486    xfac = flv2*parmDict[phfx+'Extinction']
2487    exb = 1.0
2488    if xfac > -1.:
2489        exb = 1./np.sqrt(1.+xfac)
2490    exl = 1.0
2491    if 0 < xfac <= 1.:
2492        xn = np.array([xfac**(i+1) for i in range(6)])
2493        exl += np.sum(xn*coef)
2494    elif xfac > 1.:
2495        xfac2 = 1./np.sqrt(xfac)
2496        exl = pi2*(1.-0.125/xfac)*xfac2
2497    return exb*sth2+exl*(1.-sth2)
2498   
2499def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2500    'Needs a doc string'
2501    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2502    pi2 = np.sqrt(2./np.pi)
2503    if 'T' in calcControls[hfx+'histType']:
2504        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2505        wave = refl[14+im]
2506    else:   #'C'W
2507        sth2 = sind(refl[5+im]/2.)**2
2508        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2509    c2th = 1.-2.0*sth2
2510    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2511    if 'X' in calcControls[hfx+'histType']:
2512        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2513    xfac = flv2*parmDict[phfx+'Extinction']
2514    dbde = -500.*flv2
2515    if xfac > -1.:
2516        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
2517    dlde = 0.
2518    if 0 < xfac <= 1.:
2519        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
2520        dlde = np.sum(xn*coef)
2521    elif xfac > 1.:
2522        xfac2 = 1./np.sqrt(xfac)
2523        dlde = flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
2524       
2525    return dbde*sth2+dlde*(1.-sth2)
2526   
2527def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2528    'Needs a doc string'    #need powder extinction!
2529    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
2530    if 'X' in parmDict[hfx+'Type']:
2531        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
2532    POcorr = 1.0
2533    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
2534        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2535    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
2536        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
2537    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
2538        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2539    Icorr *= POcorr
2540    AbsCorr = 1.0
2541    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
2542    Icorr *= AbsCorr
2543    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
2544    Icorr *= ExtCorr
2545    return Icorr,POcorr,AbsCorr,ExtCorr
2546   
2547def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2548    'Needs a doc string'    #need powder extinction derivs!
2549    dIdsh = 1./parmDict[hfx+'Scale']
2550    dIdsp = 1./parmDict[phfx+'Scale']
2551    if 'X' in parmDict[hfx+'Type']:
2552        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
2553        dIdPola /= pola
2554    else:       #'N'
2555        dIdPola = 0.0
2556    dFdODF = {}
2557    dFdSA = [0,0,0]
2558    dIdPO = {}
2559    if pfx+'SHorder' in parmDict:
2560        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2561        for iSH in dFdODF:
2562            dFdODF[iSH] /= odfCor
2563        for i in range(3):
2564            dFdSA[i] /= odfCor
2565    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
2566        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
2567        for iPO in dIdPO:
2568            dIdPO[iPO] /= POcorr
2569    if 'T' in parmDict[hfx+'Type']:
2570        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
2571        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
2572    else:
2573        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
2574        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
2575    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
2576       
2577def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2578    'Needs a doc string'
2579    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
2580        costh = cosd(refl[5+im]/2.)
2581        #crystallite size
2582        if calcControls[phfx+'SizeType'] == 'isotropic':
2583            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2584        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2585            H = np.array(refl[:3])
2586            P = np.array(calcControls[phfx+'SizeAxis'])
2587            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2588            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
2589            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2590        else:           #ellipsoidal crystallites
2591            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2592            H = np.array(refl[:3])
2593            lenR = G2pwd.ellipseSize(H,Sij,GB)
2594            Sgam = 1.8*wave/(np.pi*costh*lenR)
2595        #microstrain               
2596        if calcControls[phfx+'MustrainType'] == 'isotropic':
2597            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2598        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2599            H = np.array(refl[:3])
2600            P = np.array(calcControls[phfx+'MustrainAxis'])
2601            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2602            Si = parmDict[phfx+'Mustrain;i']
2603            Sa = parmDict[phfx+'Mustrain;a']
2604            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2605        else:       #generalized - P.W. Stephens model
2606            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2607            Sum = 0
2608            for i,strm in enumerate(Strms):
2609                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2610            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
2611    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
2612        #crystallite size
2613        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2614            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2615        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2616            H = np.array(refl[:3])
2617            P = np.array(calcControls[phfx+'SizeAxis'])
2618            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2619            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
2620            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2621        else:           #ellipsoidal crystallites   #OK
2622            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2623            H = np.array(refl[:3])
2624            lenR = G2pwd.ellipseSize(H,Sij,GB)
2625            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
2626        #microstrain               
2627        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
2628            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2629        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
2630            H = np.array(refl[:3])
2631            P = np.array(calcControls[phfx+'MustrainAxis'])
2632            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2633            Si = parmDict[phfx+'Mustrain;i']
2634            Sa = parmDict[phfx+'Mustrain;a']
2635            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2636        else:       #generalized - P.W. Stephens model  OK
2637            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2638            Sum = 0
2639            for i,strm in enumerate(Strms):
2640                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2641            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
2642           
2643    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
2644    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
2645    sig /= ateln2
2646    return sig,gam
2647       
2648def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2649    'Needs a doc string'
2650    gamDict = {}
2651    sigDict = {}
2652    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
2653        costh = cosd(refl[5+im]/2.)
2654        tanth = tand(refl[5+im]/2.)
2655        #crystallite size derivatives
2656        if calcControls[phfx+'SizeType'] == 'isotropic':
2657            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2658            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
2659            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
2660        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2661            H = np.array(refl[:3])
2662            P = np.array(calcControls[phfx+'SizeAxis'])
2663            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2664            Si = parmDict[phfx+'Size;i']
2665            Sa = parmDict[phfx+'Size;a']
2666            gami = 1.8*wave/(costh*np.pi*Si*Sa)
2667            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2668            Sgam = gami*sqtrm
2669            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2670            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2671            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2672            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2673            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2674            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2675        else:           #ellipsoidal crystallites
2676            const = 1.8*wave/(np.pi*costh)
2677            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2678            H = np.array(refl[:3])
2679            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2680            Sgam = const/lenR
2681            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2682                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2683                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2684        gamDict[phfx+'Size;mx'] = Sgam
2685        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
2686               
2687        #microstrain derivatives               
2688        if calcControls[phfx+'MustrainType'] == 'isotropic':
2689            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2690            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
2691            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
2692        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2693            H = np.array(refl[:3])
2694            P = np.array(calcControls[phfx+'MustrainAxis'])
2695            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2696            Si = parmDict[phfx+'Mustrain;i']
2697            Sa = parmDict[phfx+'Mustrain;a']
2698            gami = 0.018*Si*Sa*tanth/np.pi
2699            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2700            Mgam = gami/sqtrm
2701            dsi = -gami*Si*cosP**2/sqtrm**3
2702            dsa = -gami*Sa*sinP**2/sqtrm**3
2703            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2704            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2705            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2706            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2707        else:       #generalized - P.W. Stephens model
2708            const = 0.018*refl[4+im]**2*tanth/np.pi
2709            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2710            Sum = 0
2711            for i,strm in enumerate(Strms):
2712                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2713                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2714                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2715            Mgam = const*np.sqrt(Sum)
2716            for i in range(len(Strms)):
2717                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2718                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2
2719        gamDict[phfx+'Mustrain;mx'] = Mgam
2720        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2721    else:   #'T'OF - All checked & OK
2722        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2723            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2724            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
2725            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
2726        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2727            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2728            H = np.array(refl[:3])
2729            P = np.array(calcControls[phfx+'SizeAxis'])
2730            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2731            Si = parmDict[phfx+'Size;i']
2732            Sa = parmDict[phfx+'Size;a']
2733            gami = const/(Si*Sa)
2734            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2735            Sgam = gami*sqtrm
2736            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2737            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2738            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2739            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2740            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2741            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2742        else:           #OK  ellipsoidal crystallites
2743            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2744            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2745            H = np.array(refl[:3])
2746            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2747            Sgam = const/lenR
2748            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2749                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2750                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2751        gamDict[phfx+'Size;mx'] = Sgam  #OK
2752        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
2753               
2754        #microstrain derivatives               
2755        if calcControls[phfx+'MustrainType'] == 'isotropic':
2756            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2757            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
2758            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
2759        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2760            H = np.array(refl[:3])
2761            P = np.array(calcControls[phfx+'MustrainAxis'])
2762            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2763            Si = parmDict[phfx+'Mustrain;i']
2764            Sa = parmDict[phfx+'Mustrain;a']
2765            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
2766            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2767            Mgam = gami/sqtrm
2768            dsi = -gami*Si*cosP**2/sqtrm**3
2769            dsa = -gami*Sa*sinP**2/sqtrm**3
2770            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2771            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2772            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2773            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2774        else:       #generalized - P.W. Stephens model OK
2775            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2776            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
2777            Sum = 0
2778            for i,strm in enumerate(Strms):
2779                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2780                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2781                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2782            Mgam = const*np.sqrt(Sum)
2783            for i in range(len(Strms)):
2784                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2785                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2       
2786        gamDict[phfx+'Mustrain;mx'] = Mgam
2787        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2788       
2789    return sigDict,gamDict
2790       
2791def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2792    'Needs a doc string'
2793    if im:
2794        h,k,l,m = refl[:4]
2795        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2796        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
2797    else:
2798        h,k,l = refl[:3]
2799        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
2800    refl[4+im] = d
2801    if 'C' in calcControls[hfx+'histType']:
2802        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2803        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2804        if 'Bragg' in calcControls[hfx+'instType']:
2805            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2806                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2807        else:               #Debye-Scherrer - simple but maybe not right
2808            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2809    elif 'T' in calcControls[hfx+'histType']:
2810        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
2811        #do I need sample position effects - maybe?
2812    return pos
2813
2814def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2815    'Needs a doc string'
2816    dpr = 180./np.pi
2817    if im:
2818        h,k,l,m = refl[:4]
2819        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2820        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
2821        h,k,l = [h+m*vec[0],k+m*vec[1],l+m*vec[2]]          #do proj of hklm to hkl so dPdA & dPdV come out right
2822    else:
2823        m = 0
2824        h,k,l = refl[:3]       
2825        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2826    dst = np.sqrt(dstsq)
2827    dsp = 1./dst
2828    if 'C' in calcControls[hfx+'histType']:
2829        pos = refl[5+im]-parmDict[hfx+'Zero']
2830        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2831        dpdw = const*dst
2832        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
2833        dpdZ = 1.0
2834        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2835            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
2836        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2837        if 'Bragg' in calcControls[hfx+'instType']:
2838            dpdSh = -4.*shft*cosd(pos/2.0)
2839            dpdTr = -shft*sind(pos)*100.0
2840            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
2841        else:               #Debye-Scherrer - simple but maybe not right
2842            dpdXd = -shft*cosd(pos)
2843            dpdYd = -shft*sind(pos)
2844            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
2845    elif 'T' in calcControls[hfx+'histType']:
2846        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
2847        dpdZ = 1.0
2848        dpdDC = dsp
2849        dpdDA = dsp**2
2850        dpdDB = 1./dsp
2851        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2852            2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2.
2853        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2854           
2855def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2856    'Needs a doc string'
2857    laue = SGData['SGLaue']
2858    uniq = SGData['SGUniq']
2859    h,k,l = refl[:3]
2860    if laue in ['m3','m3m']:
2861        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2862            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2863    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2864        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2865    elif laue in ['3R','3mR']:
2866        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2867    elif laue in ['4/m','4/mmm']:
2868        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2869    elif laue in ['mmm']:
2870        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2871    elif laue in ['2/m']:
2872        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2873        if uniq == 'a':
2874            Dij += parmDict[phfx+'D23']*k*l
2875        elif uniq == 'b':
2876            Dij += parmDict[phfx+'D13']*h*l
2877        elif uniq == 'c':
2878            Dij += parmDict[phfx+'D12']*h*k
2879    else:
2880        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2881            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2882    if 'C' in calcControls[hfx+'histType']:
2883        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2884    else:
2885        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
2886           
2887def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2888    'Needs a doc string'
2889    laue = SGData['SGLaue']
2890    uniq = SGData['SGUniq']
2891    h,k,l = refl[:3]
2892    if laue in ['m3','m3m']:
2893        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2894            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2895    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2896        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2897    elif laue in ['3R','3mR']:
2898        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2899    elif laue in ['4/m','4/mmm']:
2900        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2901    elif laue in ['mmm']:
2902        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2903    elif laue in ['2/m']:
2904        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2905        if uniq == 'a':
2906            dDijDict[phfx+'D23'] = k*l
2907        elif uniq == 'b':
2908            dDijDict[phfx+'D13'] = h*l
2909        elif uniq == 'c':
2910            dDijDict[phfx+'D12'] = h*k
2911    else:
2912        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2913            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2914    if 'C' in calcControls[hfx+'histType']:
2915        for item in dDijDict:
2916            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2917    else:
2918        for item in dDijDict:
2919            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
2920    return dDijDict
2921   
2922def GetDij(phfx,SGData,parmDict):
2923    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
2924    return G2spc.HStrainVals(HSvals,SGData)
2925               
2926def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2927    'Compute the observed structure factors for Powder histograms'
2928    #if GSASIIpath.GetConfigValue('debug'):
2929    #    starttime = time.time(); print 'start GetFobsSq'
2930    G2mp.InitMP()
2931    histoList = Histograms.keys()
2932    histoList.sort()
2933    Ka2 = shl = lamRatio = kRatio = None
2934    for histogram in histoList:
2935        if 'PWDR' in histogram[:4]:
2936            Histogram = Histograms[histogram]
2937            hId = Histogram['hId']
2938            hfx = ':%d:'%(hId)
2939            Limits = calcControls[hfx+'Limits']
2940            if 'C' in calcControls[hfx+'histType']:
2941                shl = max(parmDict[hfx+'SH/L'],0.0005)
2942                Ka2 = False
2943                kRatio = 0.0
2944                if hfx+'Lam1' in parmDict.keys():
2945                    Ka2 = True
2946                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2947                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
2948            x,y,w,yc,yb,yd = Histogram['Data']
2949            xMask = ma.getmaskarray(x)
2950            xB = np.searchsorted(x,Limits[0])
2951            xF = np.searchsorted(x,Limits[1])
2952            ymb = np.array(y-yb)
2953            ymb = np.where(ymb,ymb,1.0)
2954            ycmb = np.array(yc-yb)
2955            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
2956            refLists = Histogram['Reflection Lists']
2957            for phase in refLists:
2958                if phase not in Phases:     #skips deleted or renamed phases silently!
2959                    continue
2960                Phase = Phases[phase]
2961                im = 0
2962                if Phase['General'].get('Modulated',False):
2963                    im = 1
2964                pId = Phase['pId']
2965                phfx = '%d:%d:'%(pId,hId)
2966                refDict = refLists[phase]
2967                sumInt = 0.0
2968                nExcl = 0
2969                useMP = G2mp.useMP and len(refDict['RefList']) > 100
2970                if useMP: # yes, create a set of initialized Python processes
2971                    MPpool = mp.Pool(G2mp.ncores,G2mp.InitFobsSqGlobals,
2972                                    [x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2])
2973                    profArgs = [[] for i in range(G2mp.ncores)]
2974                else:
2975                    G2mp.InitFobsSqGlobals(x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2)
2976                if 'C' in calcControls[hfx+'histType']:
2977                    # are we multiprocessing?
2978                    for iref,refl in enumerate(refDict['RefList']):
2979                        if useMP: 
2980                            profArgs[iref%G2mp.ncores].append((refl,iref))
2981                        else:
2982                            icod= G2mp.ComputeFobsSqCW(refl,iref)
2983                            if type(icod) is tuple:
2984                                refl[8+im] = icod[0]
2985                                sumInt += icod[1]
2986                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
2987                            elif icod == -1:
2988                                refl[3+im] *= -1
2989                                nExcl += 1
2990                            elif icod == -2:
2991                                break
2992                    if useMP:
2993                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):
2994                            sumInt += sInt
2995                            for refl8im,irefl in resList:
2996                                if refl8im is None:
2997                                    refDict['RefList'][irefl][3+im] *= -1
2998                                    nExcl += 1
2999                                else:
3000                                    refDict['RefList'][irefl][8+im] = refl8im
3001                                    if parmDict[phfx+'LeBail']:
3002                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3003                        MPpool.terminate()
3004                elif 'T' in calcControls[hfx+'histType']:
3005                    for iref,refl in enumerate(refDict['RefList']):
3006                        if useMP: 
3007                            profArgs[iref%G2mp.ncores].append((refl,iref))
3008                        else:
3009                            icod= G2mp.ComputeFobsSqTOF(refl,iref)
3010                            if type(icod) is tuple:
3011                                refl[8+im] = icod[0]
3012                                sumInt += icod[1]
3013                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3014                            elif icod == -1:
3015                                refl[3+im] *= -1
3016                                nExcl += 1
3017                            elif icod == -2:
3018                                break
3019                    if useMP:
3020                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):
3021                            sumInt += sInt
3022                            for refl8im,irefl in resList:
3023                                if refl8im is None:
3024                                    refDict['RefList'][irefl][3+im] *= -1
3025                                    nExcl += 1
3026                                else:
3027                                    refDict['RefList'][irefl][8+im] = refl8im
3028                                    if parmDict[phfx+'LeBail']:
3029                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3030                sumFo = 0.0
3031                sumdF = 0.0
3032                sumFosq = 0.0
3033                sumdFsq = 0.0
3034                for iref,refl in enumerate(refDict['RefList']):
3035                    Fo = np.sqrt(np.abs(refl[8+im]))
3036                    Fc = np.sqrt(np.abs(refl[9]+im))
3037                    sumFo += Fo
3038                    sumFosq += refl[8+im]**2
3039                    sumdF += np.abs(Fo-Fc)
3040                    sumdFsq += (refl[8+im]-refl[9+im])**2
3041                if sumFo:
3042                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3043                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3044                else:
3045                    Histogram['Residuals'][phfx+'Rf'] = 100.
3046                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
3047                Histogram['Residuals'][phfx+'sumInt'] = sumInt
3048                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])-nExcl
3049                Histogram['Residuals']['hId'] = hId
3050        elif 'HKLF' in histogram[:4]:
3051            Histogram = Histograms[histogram]
3052            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
3053    #if GSASIIpath.GetConfigValue('debug'):
3054    #    print 'end GetFobsSq t=',time.time()-starttime
3055    #    print "Rf^2",[(i,Histogram['Residuals'][i]) for i in Histogram['Residuals'] if 'Rf^2' in i]
3056
3057def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
3058    'Computes the powder pattern for a histogram based on contributions from all used phases'
3059    #if GSASIIpath.GetConfigValue('debug'):
3060    #    starttime = time.time(); print 'start getPowderProfile'
3061    G2mp.InitMP()   
3062           
3063    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
3064        U = parmDict[hfx+'U']
3065        V = parmDict[hfx+'V']
3066        W = parmDict[hfx+'W']
3067        X = parmDict[hfx+'X']
3068        Y = parmDict[hfx+'Y']
3069        tanPos = tand(refl[5+im]/2.0)
3070        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3071        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3072        sig = max(0.001,sig)
3073        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam     #save peak gamma
3074        gam = max(0.001,gam)
3075        return sig,gam
3076               
3077    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
3078        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
3079            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']/refl[4+im]**2
3080        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2
3081        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3082        sig += Ssig
3083        gam += Sgam
3084        return sig,gam
3085       
3086    def GetReflAlpBet(refl,im,hfx,parmDict):
3087        alp = parmDict[hfx+'alpha']/refl[4+im]
3088        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
3089        return alp,bet
3090   
3091    hId = Histogram['hId']
3092    hfx = ':%d:'%(hId)
3093    bakType = calcControls[hfx+'bakType']
3094    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
3095    yc = np.zeros_like(yb)
3096    cw = np.diff(x)
3097    cw = np.append(cw,cw[-1])
3098   
3099    shl = wave = kRatio = lamRatio = None # variables that are passed and must be initialized even when not needed
3100    if 'C' in calcControls[hfx+'histType']:   
3101        shl = max(parmDict[hfx+'SH/L'],0.002)
3102        Ka2 = False
3103        if hfx+'Lam1' in parmDict.keys():
3104            wave = parmDict[hfx+'Lam1']
3105            Ka2 = True
3106            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3107            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3108        else:
3109            wave = parmDict[hfx+'Lam']
3110    for phase in Histogram['Reflection Lists']:
3111        refDict = Histogram['Reflection Lists'][phase]
3112        if phase not in Phases:     #skips deleted or renamed phases silently!
3113            continue
3114        Phase = Phases[phase]
3115        pId = Phase['pId']
3116        pfx = '%d::'%(pId)
3117        phfx = '%d:%d:'%(pId,hId)
3118        hfx = ':%d:'%(hId)
3119        SGData = Phase['General']['SGData']
3120        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3121        im = 0
3122        if Phase['General'].get('Modulated',False):
3123            SSGData = Phase['General']['SSGData']
3124            im = 1  #offset in SS reflection list
3125            #??
3126        Dij = GetDij(phfx,SGData,parmDict)
3127        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
3128        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3129        if np.any(np.diag(G)<0.) or np.any(np.isnan(A)):
3130            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
3131        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3132        Vst = np.sqrt(nl.det(G))    #V*
3133        if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
3134            if im:
3135                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
3136            else:
3137                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3138        badPeak = False
3139
3140        # are we multiprocessing?
3141        useMP = G2mp.useMP and len(refDict['RefList']) > 100
3142        if useMP: # yes, create a set of initialized Python processes
3143            MPpool = mp.Pool(G2mp.ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
3144            profArgs = [[] for i in range(G2mp.ncores)]
3145        if 'C' in calcControls[hfx+'histType']:
3146            for iref,refl in enumerate(refDict['RefList']):
3147                if im:
3148                    h,k,l,m = refl[:4]
3149                else:
3150                    h,k,l = refl[:3]
3151                Uniq = np.inner(refl[:3],SGMT)
3152                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
3153                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
3154                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3155                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3156                refl[11+im] *= Vst*Lorenz
3157                 
3158                if Phase['General'].get('doPawley'):
3159                    try:
3160                        if im:
3161                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3162                        else:
3163                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3164                        refl[9+im] = parmDict[pInd]
3165                    except KeyError:
3166#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3167                        continue
3168                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
3169                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3170                iFin = np.searchsorted(x,refl[5+im]+fmax)
3171                if not iBeg+iFin:       #peak below low limit - skip peak
3172                    continue
3173