source: trunk/GSASIIstrMath.py @ 3806

Last change on this file since 3806 was 3806, checked in by toby, 3 years ago

redo parameter lookup to provide numerical deriv step size (see G2obj.getVarStep)

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