source: trunk/GSASIIstrMath.py @ 3093

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

fix multiprocessing with constraints

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