source: trunk/GSASIIstrMath.py @ 4915

Last change on this file since 4915 was 4915, checked in by vondreele, 5 months ago

add some doc strings for profile functions in G2pwd
rename getPowderProfileDervMP as getPowderProfileDerv since it is more generally used than just multiprocessing

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