source: trunk/GSASIIstrMath.py @ 4999

Last change on this file since 4999 was 4999, checked in by vondreele, 3 months ago

change incommensurate magneti str. fctr. equations
avoid search for modulation vector when no peaks were fitted.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 227.4 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2021-07-22 19:09:36 +0000 (Thu, 22 Jul 2021) $
8# $Author: vondreele $
9# $Revision: 4999 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 4999 2021-07-22 19:09:36Z 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: 4999 $")
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    Bmat = G2lat.Gmat2AB(G)[1]
1043    Mast = twopisq*np.multiply.outer(ast,ast)
1044    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1045    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1046    Ncen = len(SGData['SGCen'])
1047    Nops = len(SGMT)*Ncen
1048    if not SGData['SGFixed']:
1049        Nops *= (1+SGData['SGInv'])
1050    MFtables = calcControls['MFtables']
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    Amat,Bmat = G2lat.Gmat2AB(G)
1505    Mast = twopisq*np.multiply.outer(ast,ast)   
1506    SGInv = SGData['SGInv']
1507    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1508    Nops = len(SGMT)*(1+SGData['SGInv'])
1509    if SGData['SGGray']:
1510        Nops *= 2
1511    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1512    GamI = np.array([ops[0][3,3] for ops in SSGData['SSGOps']])
1513    if SGData['SGInv']:
1514        GamI = np.hstack((GamI,-GamI))
1515    GamI = np.hstack([GamI for cen in SGData['SGCen']])
1516    if SGData['SGGray']:
1517        GamI = np.hstack((GamI,GamI))
1518    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1519    SSCen = SSGData['SSGCen']
1520    FFtables = calcControls['FFtables']
1521    BLtables = calcControls['BLtables']
1522    MFtables = calcControls['MFtables']
1523    Flack = 1.0
1524    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1525        Flack = 1.-2.*parmDict[phfx+'Flack']
1526    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1527        GetAtomFXU(pfx,calcControls,parmDict)
1528    if not Xdata.size:          #no atoms in phase!
1529        return
1530    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1531    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)      #NB: Mmod is ReIm,Mxyz,Ntau,Natm
1532    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1533
1534    if parmDict[pfx+'isMag']:       #This part correct for making modulated mag moments on equiv atoms - Mmod matched drawing & Bilbao drawings
1535   
1536        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
1537        if SGData['SGGray']:
1538            mXYZ = np.hstack((mXYZ,mXYZ))
1539
1540        MmodAR,MmodBR,MmodAI,MmodBI = G2mth.MagMod(glTau,mXYZ,modQ,MSSdata,SGData,SSGData)  #Ntau,Nops,Natm,Mxyz cos,sin parts sum matches drawing
1541        MmodA = MmodAR+MmodBR
1542        MmodB = MmodAI+MmodBI
1543       
1544        if not SGData['SGGray']:    #for fixed Mx,My,Mz
1545            GSdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2))  #apply sym. ops.--> Natm,Nops,Nxyz
1546            if SGData['SGInv'] and not SGData['SGFixed']:   #inversion if any
1547                GSdata = np.hstack((GSdata,-GSdata))     
1548            GSdata = np.hstack([GSdata for cen in SSCen])        #dup over cell centering - Natm,Nops,Mxyz
1549            GSdata = SGData['MagMom'][nxs,:,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
1550            GSdata = np.swapaxes(GSdata,0,1)    #Nop,Natm,Mxyz
1551       
1552    FF = np.zeros(len(Tdata))
1553    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
1554        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1555    elif 'X' in calcControls[hfx+'histType']:
1556        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1557        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1558    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1559    bij = Mast*Uij
1560    blkSize = 48       #no. of reflections in a block
1561    nRef = refDict['RefList'].shape[0]
1562    SQ = 1./(2.*refDict['RefList'].T[5])**2
1563    if 'N' in calcControls[hfx+'histType']:
1564        dat = G2el.getBLvalues(BLtables)
1565        refDict['FF']['El'] = list(dat.keys())
1566        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1567        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1568        for iel,El in enumerate(refDict['FF']['El']):
1569            if El in MFtables:
1570                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1571    else:
1572        dat = G2el.getFFvalues(FFtables,0.)
1573        refDict['FF']['El'] = list(dat.keys())
1574        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1575        for iel,El in enumerate(refDict['FF']['El']):
1576            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1577#reflection processing begins here - big arrays!
1578    iBeg = 0
1579    while iBeg < nRef:
1580        iFin = min(iBeg+blkSize,nRef)
1581        mRef= iFin-iBeg
1582        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1583        H = refl.T[:4]                          #array(blkSize,4)
1584        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1585        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1586        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1587        Uniq = np.inner(H.T,SSGMT)
1588        UniqP = np.inner(HP.T,SGMT)
1589        Phi = np.inner(H.T,SSGT)
1590        if SGInv and not SGData['SGFixed']:   #if centro - expand HKL sets
1591            Uniq = np.hstack((Uniq,-Uniq))
1592            Phi = np.hstack((Phi,-Phi))
1593            UniqP = np.hstack((UniqP,-UniqP))
1594        if 'T' in calcControls[hfx+'histType']:
1595            if 'P' in calcControls[hfx+'histType']:
1596                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1597            else:
1598                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1599            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
1600            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
1601        Bab = 0.
1602        if phfx+'BabA' in parmDict:
1603            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
1604        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1605        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
1606        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1607        sinp = np.sin(phase)
1608        cosp = np.cos(phase)
1609        biso = -SQfactor*Uisodata[:,nxs]
1610        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T
1611        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1612        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1613        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1614
1615        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:           
1616           
1617            phasem = twopi*np.inner(mXYZ,HP.T).T    #2pi(Q.r)
1618            cosm = np.cos(phasem)                   #Nref,nops,natm
1619            sinm = np.sin(phasem)
1620            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1621            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Mdata*Fdata*MF/(2*Nops)     #Nref,Natm
1622            HM = np.inner(Bmat,HP.T)                #put into cartesian space X||H,Z||H*L;
1623            eM = (HM*refl.T[5]).T                   # normalize HP by d*    Nref,hkl=Unit vectors || Q
1624
1625            if not SGData['SGGray']:     #correct -fixed Mx,My,Mz contribution             
1626                fam0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*cosm[:,:,:,nxs]    #Nref,Nops,Natm,Mxyz
1627                fbm0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*sinm[:,:,:,nxs]
1628#  calc mag. structure factors; Nref,Ntau,Nops,Natm,Mxyz                           
1629            fams = TMcorr[:,nxs,nxs,:,nxs]*SGData['MagMom'][nxs,nxs,:,nxs,nxs]*np.array([np.where(H[3,i]!=0,(
1630                (MmodAR+H[3,i]*MmodBR)*cosm[i,nxs,:,:,nxs]+   
1631                GamI[nxs,:,nxs,nxs]*(-MmodAI+H[3,i]*MmodBI)*sinm[i,nxs,:,:,nxs]),
1632                0.) for i in range(mRef)])/2.          #Nref,Ntau,Nops,Natm,Mxyz
1633                       
1634            fbms = TMcorr[:,nxs,nxs,:,nxs]*SGData['MagMom'][nxs,nxs,:,nxs,nxs]*np.array([np.where(H[3,i]!=0,(
1635                (MmodAR+H[3,i]*MmodBR)*sinm[i,nxs,:,:,nxs]+   
1636                GamI[nxs,:,nxs,nxs]*(MmodAI-H[3,i]*MmodBI)*cosm[i,nxs,:,:,nxs]),
1637                0.) for i in range(mRef)])/2.          #Nref,Ntau,Nops,Natm,Mxyz
1638           
1639            if not SGData['SGGray']:
1640                fams += fam0[:,nxs,:,:,:]
1641                fbms += fbm0[:,nxs,:,:,:]
1642                               
1643#sum ops & atms                               
1644            fasm = np.sum(np.sum(fams,axis=-2),axis=-2)    #Nref,Mxyz; sum ops & atoms
1645            fbsm = np.sum(np.sum(fbms,axis=-2),axis=-2)
1646# #put into cartesian space
1647            facm = np.inner(fasm,uAmat.T)       #.T better than not
1648            fbcm = np.inner(fbsm,uAmat.T)
1649#form e.F dot product
1650            eDotFa = np.sum(eM[:,nxs,:]*facm,axis=-1)    #Nref,Ntau       
1651            eDotFb = np.sum(eM[:,nxs,:]*fbcm,axis=-1)
1652#intensity Halpern & Johnson
1653            fass = np.sum((facm-eM[:,nxs,:]*eDotFa[:,:,nxs])**2,axis=-1)
1654            fbss = np.sum((fbcm-eM[:,nxs,:]*eDotFb[:,:,nxs])**2,axis=-1)
1655## #do integration           
1656            fas = np.sum(fass*glWt[nxs,:],axis=1)
1657            fbs = np.sum(fbss*glWt[nxs,:],axis=1)
1658           
1659            refl.T[10] = fas+fbs   #Sum(fams**2,Mxyz) Re + Im
1660            refl.T[11] = atan2d(fbs,fas)
1661
1662        else:
1663            GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1664            if 'T' in calcControls[hfx+'histType']:
1665                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1666                fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1667            else:
1668                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1669                fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1670            fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1671            fbg = fb*GfpuA[0]+fa*GfpuA[1]
1672            fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1673            fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1674           
1675            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1676            refl.T[11] = atan2d(fbs[0],fas[0])  #use only tau=0; ignore f' & f"
1677        if 'P' not in calcControls[hfx+'histType']:
1678            refl.T[8] = np.copy(refl.T[10])               
1679        iBeg += blkSize
1680    return copy.deepcopy(refDict['RefList'])
1681
1682def SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1683    '''
1684    Compute super structure factors for all h,k,l,m for phase - twins only
1685    puts the result, F^2, in each ref[8+im] in refList
1686    works on blocks of 32 reflections for speed
1687    input:
1688   
1689    :param dict refDict: where
1690        'RefList' list where each ref = h,k,l,m,it,d,...
1691        'FF' dict of form factors - filed in below
1692    :param np.array G:      reciprocal metric tensor
1693    :param str pfx:    phase id string
1694    :param dict SGData: space group info. dictionary output from SpcGroup
1695    :param dict calcControls:
1696    :param dict ParmDict:
1697
1698    '''
1699    phfx = pfx.split(':')[0]+hfx
1700    ast = np.sqrt(np.diag(G))
1701    Mast = twopisq*np.multiply.outer(ast,ast)   
1702    SGInv = SGData['SGInv']
1703    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1704    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1705    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1706    FFtables = calcControls['FFtables']
1707    BLtables = calcControls['BLtables']
1708    MFtables = calcControls['MFtables']
1709    Flack = 1.0
1710    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1711        Flack = 1.-2.*parmDict[phfx+'Flack']
1712    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])    #4D?
1713    TwDict = refDict.get('TwDict',{})           
1714    if 'S' in calcControls[hfx+'histType']:
1715        NTL = calcControls[phfx+'NTL']
1716        NM = calcControls[phfx+'TwinNMN']+1
1717        TwinLaw = calcControls[phfx+'TwinLaw']  #this'll have to be 4D also...
1718        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1719        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1720    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1721        GetAtomFXU(pfx,calcControls,parmDict)
1722    if not Xdata.size:          #no atoms in phase!
1723        return
1724    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1725    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1726    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1727    FF = np.zeros(len(Tdata))
1728    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
1729        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1730    elif 'X' in calcControls[hfx+'histType']:
1731        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1732        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1733    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1734    bij = Mast*Uij
1735    blkSize = 32       #no. of reflections in a block
1736    nRef = refDict['RefList'].shape[0]
1737    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
1738        SQ = 1./(2.*refDict['RefList'].T[5])**2
1739        if 'N' in calcControls[hfx+'histType']:
1740            dat = G2el.getBLvalues(BLtables)
1741            refDict['FF']['El'] = list(dat.keys())
1742            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1743            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1744            for iel,El in enumerate(refDict['FF']['El']):
1745                if El in MFtables:
1746                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1747        else:
1748            dat = G2el.getFFvalues(FFtables,0.)
1749            refDict['FF']['El'] = list(dat.keys())
1750            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1751            for iel,El in enumerate(refDict['FF']['El']):
1752                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1753#    time0 = time.time()
1754#reflection processing begins here - big arrays!
1755    iBeg = 0
1756    while iBeg < nRef:
1757        iFin = min(iBeg+blkSize,nRef)
1758        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1759        H = refl[:,:4]                          #array(blkSize,4)
1760        H3 = refl[:,:3]
1761        HP = H[:,:3]+modQ[nxs,:]*H[:,3:]        #projected hklm to hkl
1762        HP = np.inner(HP,TwinLaw)             #array(blkSize,nTwins,4)
1763        H3 = np.inner(H3,TwinLaw)       
1764        TwMask = np.any(HP,axis=-1)
1765        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1766            for ir in range(blkSize):
1767                iref = ir+iBeg
1768                if iref in TwDict:
1769                    for i in TwDict[iref]:
1770                        for n in range(NTL):
1771                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1772                            H3[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1773            TwMask = np.any(HP,axis=-1)
1774        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1775        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1776        Uniq = np.inner(H,SSGMT)
1777        Uniq3 = np.inner(H3,SGMT)
1778        UniqP = np.inner(HP,SGMT)
1779        Phi = np.inner(H,SSGT)
1780        if SGInv:   #if centro - expand HKL sets
1781            Uniq = np.hstack((Uniq,-Uniq))
1782            Uniq3 = np.hstack((Uniq3,-Uniq3))
1783            Phi = np.hstack((Phi,-Phi))
1784            UniqP = np.hstack((UniqP,-UniqP))
1785        if 'T' in calcControls[hfx+'histType']:
1786            if 'P' in calcControls[hfx+'histType']:
1787                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1788            else:
1789                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1790            FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1791            FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1792        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw))
1793        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1794        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1795        phase = twopi*(np.inner(Uniq3,(dXdata.T+Xdata.T))-Phi[:,nxs,:,nxs])
1796        sinp = np.sin(phase)
1797        cosp = np.cos(phase)
1798        biso = -SQfactor*Uisodata[:,nxs]
1799        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1800        HbH = -np.sum(UniqP[:,:,:,nxs]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1801        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1802        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1803        if 'T' in calcControls[hfx+'histType']:
1804            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1805            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1806        else:
1807            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1808            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1809        GfpuA = G2mth.ModulationTw(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1810        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1811        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1812        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1813        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1814        refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
1815        refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
1816            np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
1817        refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
1818        iBeg += blkSize
1819#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1820
1821def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1822    '''
1823    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
1824    Only Fourier component are done analytically here
1825    input:
1826   
1827    :param dict refDict: where
1828        'RefList' list where each ref = h,k,l,m,it,d,...
1829        'FF' dict of form factors - filled in below
1830    :param int im: = 1 (could be eliminated)
1831    :param np.array G:      reciprocal metric tensor
1832    :param str hfx:    histogram id string
1833    :param str pfx:    phase id string
1834    :param dict SGData: space group info. dictionary output from SpcGroup
1835    :param dict SSGData: super space group info.
1836    :param dict calcControls:
1837    :param dict ParmDict:
1838   
1839    :returns: dict dFdvDict: dictionary of derivatives
1840    '''
1841    phfx = pfx.split(':')[0]+hfx
1842    ast = np.sqrt(np.diag(G))
1843    Mast = twopisq*np.multiply.outer(ast,ast)
1844    SGInv = SGData['SGInv']
1845    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1846    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1847    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1848    FFtables = calcControls['FFtables']
1849    BLtables = calcControls['BLtables']
1850    nRef = len(refDict['RefList'])
1851    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1852        GetAtomFXU(pfx,calcControls,parmDict)
1853    if not Xdata.size:          #no atoms in phase!
1854        return {}
1855    mSize = len(Mdata)  #no. atoms
1856    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1857    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
1858    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
1859    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1860    FF = np.zeros(len(Tdata))
1861    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
1862        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1863    elif 'X' in calcControls[hfx+'histType']:
1864        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1865        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1866    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1867    bij = Mast*Uij
1868    if not len(refDict['FF']):
1869        if 'N' in calcControls[hfx+'histType']:
1870            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1871        else:
1872            dat = G2el.getFFvalues(FFtables,0.)       
1873        refDict['FF']['El'] = list(dat.keys())
1874        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1875    dFdvDict = {}
1876    dFdfr = np.zeros((nRef,mSize))
1877    dFdx = np.zeros((nRef,mSize,3))
1878    dFdui = np.zeros((nRef,mSize))
1879    dFdua = np.zeros((nRef,mSize,6))
1880    dFdbab = np.zeros((nRef,2))
1881    dFdfl = np.zeros((nRef))
1882    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1883    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1884    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1885    Flack = 1.0
1886    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1887        Flack = 1.-2.*parmDict[phfx+'Flack']
1888    time0 = time.time()
1889    nRef = len(refDict['RefList'])/100
1890    for iref,refl in enumerate(refDict['RefList']):
1891        if 'T' in calcControls[hfx+'histType']:
1892            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1893        H = np.array(refl[:4])
1894        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
1895        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1896        SQfactor = 8.0*SQ*np.pi**2
1897        Bab = 0.0
1898        if phfx+'BabA' in parmDict:
1899            dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1900            Bab = parmDict[phfx+'BabA']*dBabdA
1901        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1902        FF = refDict['FF']['FF'][iref].T[Tindx]
1903        Uniq = np.inner(H,SSGMT)
1904        Phi = np.inner(H,SSGT)
1905        UniqP = np.inner(HP,SGMT)
1906        if SGInv:   #if centro - expand HKL sets
1907            Uniq = np.vstack((Uniq,-Uniq))
1908            Phi = np.hstack((Phi,-Phi))
1909            UniqP = np.vstack((UniqP,-UniqP))
1910        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1911        sinp = np.sin(phase)
1912        cosp = np.cos(phase)
1913        occ = Mdata*Fdata/Uniq.shape[0]
1914        biso = -SQfactor*Uisodata[:,nxs]
1915        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0],axis=1).T    #ops x atoms
1916        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
1917        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1918        Hij = np.array([G2lat.UijtoU6(uij) for uij in Hij])                     #atoms x 6
1919        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1920        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1921        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
1922        fotp = FPP*Tcorr            #ops x atoms
1923        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1924        dGdf,dGdx,dGdu = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1925        # GfpuA is 2 x ops x atoms
1926        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1927        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
1928        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1929        fag = fa*GfpuA[0]-fb*GfpuA[1]
1930        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1931       
1932        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
1933        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1934        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
1935        fbx = np.array([fot*cosp,-fotp*sinp])
1936        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1937        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1938        #sum below is over Uniq
1939        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
1940        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
1941        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1942        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1943        dfadui = np.sum(-SQfactor*fag,axis=1)
1944        dfbdui = np.sum(-SQfactor*fbg,axis=1)
1945        dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
1946        dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
1947        dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
1948        dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
1949        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1950        dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1951        dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1952        dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1953        dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1954        dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
1955        dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
1956        if not SGData['SGInv']:   #Flack derivative
1957            dfadfl = np.sum(-FPP*Tcorr*sinp)
1958            dfbdfl = np.sum(FPP*Tcorr*cosp)
1959        else:
1960            dfadfl = 1.0
1961            dfbdfl = 1.0
1962        SA = fas[0]+fas[1]      #float = A+A'
1963        SB = fbs[0]+fbs[1]      #float = B+B'
1964        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
1965            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1966            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
1967                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1968            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
1969                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1970            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
1971                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1972            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
1973                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1974            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
1975                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
1976            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
1977                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
1978            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
1979                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
1980        else:                       #OK, I think
1981            dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
1982            dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
1983            dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
1984            dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
1985            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1986                           
1987            dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
1988            dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
1989            dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
1990        if phfx+'BabA' in parmDict:
1991            dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1992                2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1993        #loop over atoms - each dict entry is list of derivatives for all the reflections
1994        if not iref%100 :
1995            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
1996    for i in range(len(Mdata)):     #loop over atoms
1997        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1998        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1999        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2000        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2001        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2002        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2003        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2004        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2005        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
2006        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
2007        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
2008        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2009            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2010            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2011        nx = 0
2012        if waveTypes[i] in ['Block','ZigZag']:
2013            nx = 1 
2014        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2015            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2016            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2017            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2018            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2019            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2020            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2021        for j in range(USSdata.shape[1]):       #loop over waves
2022            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2023            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2024            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2025            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
2026            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
2027            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
2028            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2029            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2030            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2031            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
2032            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
2033            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
2034           
2035    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
2036    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2037    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2038    return dFdvDict
2039
2040def SStructureFactorDerv2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2041    '''
2042    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
2043    input:
2044   
2045    :param dict refDict: where
2046        'RefList' list where each ref = h,k,l,m,it,d,...
2047        'FF' dict of form factors - filled in below
2048    :param int im: = 1 (could be eliminated)
2049    :param np.array G:      reciprocal metric tensor
2050    :param str hfx:    histogram id string
2051    :param str pfx:    phase id string
2052    :param dict SGData: space group info. dictionary output from SpcGroup
2053    :param dict SSGData: super space group info.
2054    :param dict calcControls:
2055    :param dict ParmDict:
2056   
2057    :returns: dict dFdvDict: dictionary of derivatives
2058    '''
2059
2060    trefDict = copy.deepcopy(refDict)
2061    dM = 1.e-4
2062    dFdvDict = {}
2063    for parm in parmDict:
2064        if ':' not in parm:
2065            continue
2066        if parm.split(':')[2] in ['Tmin','Tmax','Xmax','Ymax','Zmax','Fzero','Fwid',
2067            'MXsin','MXcos','MYsin','MYcos','MZsin','MZcos','AMx','AMy','AMz',]:
2068            parmDict[parm] += dM
2069            prefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2070            parmDict[parm] -= 2*dM
2071            mrefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2072            parmDict[parm] += dM
2073            dFdvDict[parm] = (prefList[:,9+im]-mrefList[:,9+im])/(2.*dM)
2074    return dFdvDict
2075   
2076def SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2077    'Needs a doc string'
2078    phfx = pfx.split(':')[0]+hfx
2079    ast = np.sqrt(np.diag(G))
2080    Mast = twopisq*np.multiply.outer(ast,ast)
2081    SGInv = SGData['SGInv']
2082    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2083    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2084    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2085    FFtables = calcControls['FFtables']
2086    BLtables = calcControls['BLtables']
2087    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
2088    TwDict = refDict.get('TwDict',{})           
2089    if 'S' in calcControls[hfx+'histType']:
2090        NTL = calcControls[phfx+'NTL']
2091        NM = calcControls[phfx+'TwinNMN']+1
2092        TwinLaw = calcControls[phfx+'TwinLaw']
2093        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
2094    nTwin = len(TwinLaw)       
2095    nRef = len(refDict['RefList'])
2096    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
2097        GetAtomFXU(pfx,calcControls,parmDict)
2098    if not Xdata.size:          #no atoms in phase!
2099        return {}
2100    mSize = len(Mdata)  #no. atoms
2101    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
2102    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)     #NB: Mmod is ReIm,Mxyz,Ntau,Natm
2103    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
2104    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2105    FF = np.zeros(len(Tdata))
2106    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
2107        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
2108    elif 'X' in calcControls[hfx+'histType']:
2109        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2110        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2111    Uij = np.array(G2lat.U6toUij(Uijdata)).T
2112    bij = Mast*Uij
2113    if not len(refDict['FF']):
2114        if 'N' in calcControls[hfx+'histType']:
2115            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
2116        else:
2117            dat = G2el.getFFvalues(FFtables,0.)       
2118        refDict['FF']['El'] = list(dat.keys())
2119        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
2120    dFdvDict = {}
2121    dFdfr = np.zeros((nRef,nTwin,mSize))
2122    dFdx = np.zeros((nRef,nTwin,mSize,3))
2123    dFdui = np.zeros((nRef,nTwin,mSize))
2124    dFdua = np.zeros((nRef,nTwin,mSize,6))
2125    dFdbab = np.zeros((nRef,nTwin,2))
2126    dFdtw = np.zeros((nRef,nTwin))
2127    dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
2128    dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
2129    dFdGz = np.zeros((nRef,nTwin,mSize,5))
2130    dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
2131    Flack = 1.0
2132    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
2133        Flack = 1.-2.*parmDict[phfx+'Flack']
2134    time0 = time.time()
2135    nRef = len(refDict['RefList'])/100
2136    for iref,refl in enumerate(refDict['RefList']):
2137        if 'T' in calcControls[hfx+'histType']:
2138            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
2139        H = np.array(refl[:4])
2140        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
2141        H = np.inner(H.T,TwinLaw)   #maybe array(4,nTwins) or (4)
2142        TwMask = np.any(H,axis=-1)
2143        if TwinLaw.shape[0] > 1 and TwDict:
2144            if iref in TwDict:
2145                for i in TwDict[iref]:
2146                    for n in range(NTL):
2147                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
2148            TwMask = np.any(H,axis=-1)
2149        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
2150        SQfactor = 8.0*SQ*np.pi**2
2151        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2152        Bab = parmDict[phfx+'BabA']*dBabdA
2153        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
2154        FF = refDict['FF']['FF'][iref].T[Tindx]
2155        Uniq = np.inner(H,SSGMT)
2156        Phi = np.inner(H,SSGT)
2157        UniqP = np.inner(HP,SGMT)
2158        if SGInv:   #if centro - expand HKL sets
2159            Uniq = np.vstack((Uniq,-Uniq))
2160            Phi = np.hstack((Phi,-Phi))
2161            UniqP = np.vstack((UniqP,-UniqP))
2162        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
2163        sinp = np.sin(phase)
2164        cosp = np.cos(phase)
2165        occ = Mdata*Fdata/Uniq.shape[0]
2166        biso = -SQfactor*Uisodata[:,nxs]
2167        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
2168        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
2169        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
2170        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(nTwin,-1,6)))
2171        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
2172        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
2173        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
2174        fotp = FPP*Tcorr            #ops x atoms
2175        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
2176        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
2177        # GfpuA is 2 x ops x atoms
2178        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
2179        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
2180        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
2181        fag = fa*GfpuA[0]-fb*GfpuA[1]
2182        fbg = fb*GfpuA[0]+fa*GfpuA[1]
2183       
2184        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
2185        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
2186        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x twin x ops x atoms
2187        fbx = np.array([fot*cosp,-fotp*sinp])
2188        fax = fax*GfpuA[0]-fbx*GfpuA[1]
2189        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
2190        #sum below is over Uniq
2191        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
2192        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
2193        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
2194        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
2195        dfadui = np.sum(-SQfactor*fag,axis=1)
2196        dfbdui = np.sum(-SQfactor*fbg,axis=1)
2197        dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2198        dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
2199        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2200        dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2201        # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
2202        dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2203        dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2204        dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2205        dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2206        dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2207        dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2208        dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2209        dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2210#        GSASIIpath.IPyBreak()
2211        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
2212        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
2213        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
2214        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)]
2215        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)]
2216        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)]
2217        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)]
2218        dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
2219
2220        dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
2221        dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
2222        dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
2223        dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
2224#            GSASIIpath.IPyBreak()
2225        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
2226            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2227        #loop over atoms - each dict entry is list of derivatives for all the reflections
2228        if not iref%100 :
2229            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
2230    for i in range(len(Mdata)):     #loop over atoms
2231        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2232        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2233        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2234        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2235        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2236        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2237        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2238        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2239        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
2240        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
2241        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
2242        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2243            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2244            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2245        nx = 0
2246        if waveTypes[i] in ['Block','ZigZag']:
2247            nx = 1 
2248            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
2249            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
2250            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
2251            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
2252            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
2253        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2254            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2255            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2256            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2257            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2258            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2259            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2260        for j in range(USSdata.shape[1]):       #loop over waves
2261            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2262            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2263            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2264            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
2265            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
2266            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
2267            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2268            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2269            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2270            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
2271            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
2272            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
2273           
2274#        GSASIIpath.IPyBreak()
2275    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2276    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2277    return dFdvDict
2278   
2279def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
2280    ''' Single crystal extinction function; returns extinction & derivative
2281    '''
2282    extCor = 1.0
2283    dervDict = {}
2284    dervCor = 1.0
2285    if calcControls[phfx+'EType'] != 'None':
2286        SQ = 1/(4.*ref[4+im]**2)
2287        if 'C' in parmDict[hfx+'Type']:           
2288            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
2289        else:   #'T'
2290            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
2291        if 'SXC' in parmDict[hfx+'Type']:
2292            AV = 7.9406e5/parmDict[pfx+'Vol']**2
2293            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2294            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
2295            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
2296        elif 'SNT' in parmDict[hfx+'Type']:
2297            AV = 1.e7/parmDict[pfx+'Vol']**2
2298            PL = SQ
2299            PLZ = AV*ref[9+im]*ref[12+im]**2
2300        elif 'SNC' in parmDict[hfx+'Type']:
2301            AV = 1.e7/parmDict[pfx+'Vol']**2
2302            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2303            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
2304           
2305        if 'Primary' in calcControls[phfx+'EType']:
2306            PLZ *= 1.5
2307        else:
2308            if 'C' in parmDict[hfx+'Type']:
2309                PLZ *= calcControls[phfx+'Tbar']
2310            else: #'T'
2311                PLZ *= ref[13+im]      #t-bar
2312        if 'Primary' in calcControls[phfx+'EType']:
2313            PLZ *= 1.5
2314            PSIG = parmDict[phfx+'Ep']
2315        elif 'I & II' in calcControls[phfx+'EType']:
2316            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
2317        elif 'Type II' in calcControls[phfx+'EType']:
2318            PSIG = parmDict[phfx+'Es']
2319        else:       # 'Secondary Type I'
2320            PSIG = parmDict[phfx+'Eg']/PL
2321           
2322        AG = 0.58+0.48*cos2T+0.24*cos2T**2
2323        AL = 0.025+0.285*cos2T
2324        BG = 0.02-0.025*cos2T
2325        BL = 0.15-0.2*(0.75-cos2T)**2
2326        if cos2T < 0.:
2327            BL = -0.45*cos2T
2328        CG = 2.
2329        CL = 2.
2330        PF = PLZ*PSIG
2331       
2332        if 'Gaussian' in calcControls[phfx+'EApprox']:
2333            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
2334            extCor = np.sqrt(PF4)
2335            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
2336        else:
2337            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
2338            extCor = np.sqrt(PF4)
2339            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
2340
2341        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
2342        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
2343            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
2344        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
2345            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
2346        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
2347            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
2348               
2349    return 1./extCor,dervDict,dervCor
2350   
2351def Dict2Values(parmdict, varylist):
2352    '''Use before call to leastsq to setup list of values for the parameters
2353    in parmdict, as selected by key in varylist'''
2354    return [parmdict[key] for key in varylist] 
2355   
2356def Values2Dict(parmdict, varylist, values):
2357    ''' Use after call to leastsq to update the parameter dictionary with
2358    values corresponding to keys in varylist'''
2359    parmdict.update(zip(varylist,values))
2360   
2361def GetNewCellParms(parmDict,varyList):
2362    '''Compute unit cell tensor terms from varied Aij and Dij values.
2363    Terms are included in the dict only if Aij or Dij is varied.
2364    '''
2365    newCellDict = {}
2366    Anames = ['A'+str(i) for i in range(6)]
2367    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
2368    for item in varyList:
2369        keys = item.split(':')
2370        if keys[2] in Ddict:
2371            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
2372            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
2373            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2374    return newCellDict          # is e.g. {'0::D11':A0-D11}
2375   
2376def ApplyXYZshifts(parmDict,varyList):
2377    '''
2378    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2379   
2380    :param dict parmDict: parameter dictionary
2381    :param list varyList: list of variables (not used!)
2382    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
2383
2384    '''
2385    newAtomDict = {}
2386    for item in parmDict:
2387        if 'dA' in item:
2388            parm = ''.join(item.split('d'))
2389            parmDict[parm] += parmDict[item]
2390            newAtomDict[item] = [parm,parmDict[parm]]
2391    return newAtomDict
2392   
2393def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2394    'Spherical harmonics texture'
2395    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2396    if 'T' in calcControls[hfx+'histType']:
2397        tth = parmDict[hfx+'2-theta']
2398    else:
2399        tth = refl[5+im]
2400    odfCor = 1.0
2401    H = refl[:3]
2402    cell = G2lat.Gmat2cell(g)
2403    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2404    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2405    phi,beta = G2lat.CrsAng(H,cell,SGData)
2406    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2407    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2408    for item in SHnames:
2409        L,M,N = eval(item.strip('C'))
2410        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2411        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2412        Lnorm = G2lat.Lnorm(L)
2413        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2414    return odfCor
2415   
2416def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2417    'Spherical harmonics texture derivatives'
2418    if 'T' in calcControls[hfx+'histType']:
2419        tth = parmDict[hfx+'2-theta']
2420    else:
2421        tth = refl[5+im]
2422    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2423    odfCor = 1.0
2424    dFdODF = {}
2425    dFdSA = [0,0,0]
2426    H = refl[:3]
2427    cell = G2lat.Gmat2cell(g)
2428    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2429    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2430    phi,beta = G2lat.CrsAng(H,cell,SGData)
2431    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
2432    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2433    for item in SHnames:
2434        L,M,N = eval(item.strip('C'))
2435        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2436        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2437        Lnorm = G2lat.Lnorm(L)
2438        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2439        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2440        for i in range(3):
2441            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2442    return odfCor,dFdODF,dFdSA
2443   
2444def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2445    'spherical harmonics preferred orientation (cylindrical symmetry only)'
2446    if 'T' in calcControls[hfx+'histType']:
2447        tth = parmDict[hfx+'2-theta']
2448    else:
2449        tth = refl[5+im]
2450    odfCor = 1.0
2451    H = refl[:3]
2452    cell = G2lat.Gmat2cell(g)
2453    Sangls = [0.,0.,0.]
2454    if 'Bragg' in calcControls[hfx+'instType']:
2455        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2456        IFCoup = True
2457    else:
2458        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2459        IFCoup = False
2460    phi,beta = G2lat.CrsAng(H,cell,SGData)
2461    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2462    SHnames = calcControls[phfx+'SHnames']
2463    for item in SHnames:
2464        L,N = eval(item.strip('C'))
2465        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2466        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2467        Lnorm = G2lat.Lnorm(L)
2468        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2469    return np.squeeze(odfCor)
2470   
2471def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2472    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
2473    if 'T' in calcControls[hfx+'histType']:
2474        tth = parmDict[hfx+'2-theta']
2475    else:
2476        tth = refl[5+im]
2477    odfCor = 1.0
2478    dFdODF = {}
2479    H = refl[:3]
2480    cell = G2lat.Gmat2cell(g)
2481    Sangls = [0.,0.,0.]
2482    if 'Bragg' in calcControls[hfx+'instType']:
2483        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2484        IFCoup = True
2485    else:
2486        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2487        IFCoup = False
2488    phi,beta = G2lat.CrsAng(H,cell,SGData)
2489    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2490    SHnames = calcControls[phfx+'SHnames']
2491    for item in SHnames:
2492        L,N = eval(item.strip('C'))
2493        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2494        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2495        Lnorm = G2lat.Lnorm(L)
2496        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2497        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
2498    return odfCor,dFdODF
2499   
2500def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2501    'March-Dollase preferred orientation correction'
2502    POcorr = 1.0
2503    MD = parmDict[phfx+'MD']
2504    if MD != 1.0:
2505        MDAxis = calcControls[phfx+'MDAxis']
2506        sumMD = 0
2507        for H in uniq:           
2508            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2509            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2510            sumMD += A**3
2511        POcorr = sumMD/len(uniq)
2512    return POcorr
2513   
2514def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2515    'Needs a doc string'
2516    POcorr = 1.0
2517    POderv = {}
2518    if calcControls[phfx+'poType'] == 'MD':
2519        MD = parmDict[phfx+'MD']
2520        MDAxis = calcControls[phfx+'MDAxis']
2521        sumMD = 0
2522        sumdMD = 0
2523        for H in uniq:           
2524            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2525            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2526            sumMD += A**3
2527            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2528        POcorr = sumMD/len(uniq)
2529        POderv[phfx+'MD'] = sumdMD/len(uniq)
2530    else:   #spherical harmonics
2531        if calcControls[phfx+'SHord']:
2532            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2533    return POcorr,POderv
2534   
2535def GetAbsorb(refl,im,hfx,calcControls,parmDict):
2536    'Needs a doc string'
2537    if 'Debye' in calcControls[hfx+'instType']:
2538        if 'T' in calcControls[hfx+'histType']:
2539            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2540        else:
2541            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2542    else:
2543        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
2544   
2545def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
2546    'Needs a doc string'
2547    if 'Debye' in calcControls[hfx+'instType']:
2548        if 'T' in calcControls[hfx+'histType']:
2549            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2550        else:
2551            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2552    else:
2553        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
2554       
2555def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2556    'Needs a doc string'
2557    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2558    pi2 = np.sqrt(2./np.pi)
2559    if 'T' in calcControls[hfx+'histType']:
2560        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2561        wave = refl[14+im]
2562    else:   #'C'W
2563        sth2 = sind(refl[5+im]/2.)**2
2564        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2565    c2th = 1.-2.0*sth2
2566    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2567    if 'X' in calcControls[hfx+'histType']:
2568        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2569    xfac = flv2*parmDict[phfx+'Extinction']
2570    exb = 1.0
2571    if xfac > -1.:
2572        exb = 1./np.sqrt(1.+xfac)
2573    exl = 1.0
2574    if 0 < xfac <= 1.:
2575        xn = np.array([xfac**(i+1) for i in range(6)])
2576        exl += np.sum(xn*coef)
2577    elif xfac > 1.:
2578        xfac2 = 1./np.sqrt(xfac)
2579        exl = pi2*(1.-0.125/xfac)*xfac2
2580    return exb*sth2+exl*(1.-sth2)
2581   
2582def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2583    'Needs a doc string'
2584    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2585    pi2 = np.sqrt(2./np.pi)
2586    if 'T' in calcControls[hfx+'histType']:
2587        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2588        wave = refl[14+im]
2589    else:   #'C'W
2590        sth2 = sind(refl[5+im]/2.)**2
2591        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2592    c2th = 1.-2.0*sth2
2593    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2594    if 'X' in calcControls[hfx+'histType']:
2595        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2596    xfac = flv2*parmDict[phfx+'Extinction']
2597    dbde = -500.*flv2
2598    if xfac > -1.:
2599        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
2600    dlde = 0.
2601    if 0 < xfac <= 1.:
2602        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
2603        dlde = np.sum(xn*coef)/xfac
2604    elif xfac > 1.:
2605        xfac2 = 1./np.sqrt(xfac)
2606        dlde = 0.5*flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
2607       
2608    return dbde*sth2+dlde*(1.-sth2)
2609   
2610def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2611    'Needs a doc string'    #need powder extinction!
2612    parmDict[phfx+'Scale'] = max(1.e-12,parmDict[phfx+'Scale'])                      #put floor on phase fraction scale
2613    parmDict[hfx+'Scale'] = max(1.e-12,parmDict[hfx+'Scale'])                        #put floor on histogram scale
2614    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
2615    if 'X' in parmDict[hfx+'Type']:
2616        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
2617    POcorr = 1.0
2618    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
2619        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2620    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
2621        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
2622    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
2623        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2624    Icorr *= POcorr
2625    AbsCorr = 1.0
2626    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
2627    Icorr *= AbsCorr
2628    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
2629    Icorr *= ExtCorr
2630    return Icorr,POcorr,AbsCorr,ExtCorr
2631   
2632def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2633    'Needs a doc string'    #need powder extinction derivs!
2634    dIdsh = 1./parmDict[hfx+'Scale']
2635    dIdsp = 1./parmDict[phfx+'Scale']
2636    if 'X' in parmDict[hfx+'Type']:
2637        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
2638        dIdPola /= pola
2639    else:       #'N'
2640        dIdPola = 0.0
2641    dFdODF = {}
2642    dFdSA = [0,0,0]
2643    dIdPO = {}
2644    if pfx+'SHorder' in parmDict:
2645        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2646        for iSH in dFdODF:
2647            dFdODF[iSH] /= odfCor
2648        for i in range(3):
2649            dFdSA[i] /= odfCor
2650    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
2651        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
2652        for iPO in dIdPO:
2653            dIdPO[iPO] /= POcorr
2654    if 'T' in parmDict[hfx+'Type']:
2655        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
2656        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
2657    else:
2658        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
2659        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
2660    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
2661       
2662def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2663    'Needs a doc string'
2664    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:     #All checked & OK
2665        costh = cosd(refl[5+im]/2.)
2666        #crystallite size
2667        if calcControls[phfx+'SizeType'] == 'isotropic':
2668            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2669        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2670            H = np.array(refl[:3])
2671            P = np.array(calcControls[phfx+'SizeAxis'])
2672            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2673            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
2674            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2675        else:           #ellipsoidal crystallites
2676            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2677            H = np.array(refl[:3])
2678            lenR = G2pwd.ellipseSize(H,Sij,GB)
2679            Sgam = 1.8*wave/(np.pi*costh*lenR)
2680        #microstrain               
2681        if calcControls[phfx+'MustrainType'] == 'isotropic':
2682            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2683        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2684            H = np.array(refl[:3])
2685            P = np.array(calcControls[phfx+'MustrainAxis'])
2686            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2687            Si = parmDict[phfx+'Mustrain;i']
2688            Sa = parmDict[phfx+'Mustrain;a']
2689            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2690        else:       #generalized - P.W. Stephens model
2691            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2692            Sum = 0
2693            for i,strm in enumerate(Strms):
2694                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2695            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
2696    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
2697        #crystallite size
2698        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2699            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2700        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2701            H = np.array(refl[:3])
2702            P = np.array(calcControls[phfx+'SizeAxis'])
2703            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2704            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
2705            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2706        else:           #ellipsoidal crystallites   #OK
2707            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2708            H = np.array(refl[:3])
2709            lenR = G2pwd.ellipseSize(H,Sij,GB)
2710            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
2711        #microstrain               
2712        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
2713            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2714        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
2715            H = np.array(refl[:3])
2716            P = np.array(calcControls[phfx+'MustrainAxis'])
2717            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2718            Si = parmDict[phfx+'Mustrain;i']
2719            Sa = parmDict[phfx+'Mustrain;a']
2720            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2721        else:       #generalized - P.W. Stephens model  OK
2722            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2723            Sum = 0
2724            for i,strm in enumerate(Strms):
2725                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2726            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
2727           
2728    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
2729    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
2730    sig /= ateln2
2731    return sig,gam
2732       
2733def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2734    'Needs a doc string'
2735    gamDict = {}
2736    sigDict = {}
2737    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:         #All checked & OK
2738        costh = cosd(refl[5+im]/2.)
2739        tanth = tand(refl[5+im]/2.)
2740        #crystallite size derivatives
2741        if calcControls[phfx+'SizeType'] == 'isotropic':
2742            Sgam = 1.8*wave/(np.pi*costh*parmDict[phfx+'Size;i'])
2743            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh*parmDict[phfx+'Size;i']**2)
2744            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
2745        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2746            H = np.array(refl[:3])
2747            P = np.array(calcControls[phfx+'SizeAxis'])
2748            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2749            Si = parmDict[phfx+'Size;i']
2750            Sa = parmDict[phfx+'Size;a']
2751            gami = 1.8*wave/(costh*np.pi*Si*Sa)
2752            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2753            Sgam = gami*sqtrm
2754            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2755            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2756            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2757            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2758            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2759            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2760        else:           #ellipsoidal crystallites
2761            const = 1.8*wave/(np.pi*costh)
2762            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2763            H = np.array(refl[:3])
2764            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2765            Sgam = const/lenR
2766            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2767                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2768                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2769        gamDict[phfx+'Size;mx'] = Sgam
2770        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
2771               
2772        #microstrain derivatives               
2773        if calcControls[phfx+'MustrainType'] == 'isotropic':
2774            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2775            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
2776            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
2777        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2778            H = np.array(refl[:3])
2779            P = np.array(calcControls[phfx+'MustrainAxis'])
2780            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2781            Si = parmDict[phfx+'Mustrain;i']
2782            Sa = parmDict[phfx+'Mustrain;a']
2783            gami = 0.018*Si*Sa*tanth/np.pi
2784            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2785            Mgam = gami/sqtrm
2786            dsi = -gami*Si*cosP**2/sqtrm**3
2787            dsa = -gami*Sa*sinP**2/sqtrm**3
2788            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2789            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2790            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2791            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2792        else:       #generalized - P.W. Stephens model
2793            const = 0.018*refl[4+im]**2*tanth/np.pi
2794            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2795            Sum = 0
2796            for i,strm in enumerate(Strms):
2797                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2798                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2799                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2800            Mgam = const*np.sqrt(Sum)
2801            for i in range(len(Strms)):
2802                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2803                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2
2804        gamDict[phfx+'Mustrain;mx'] = Mgam
2805        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2806    else:   #'T'OF - All checked & OK
2807        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2808            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2809            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
2810            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
2811        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2812            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2813            H = np.array(refl[:3])
2814            P = np.array(calcControls[phfx+'SizeAxis'])
2815            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2816            Si = parmDict[phfx+'Size;i']
2817            Sa = parmDict[phfx+'Size;a']
2818            gami = const/(Si*Sa)
2819            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2820            Sgam = gami*sqtrm
2821            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2822            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2823            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2824            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2825            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2826            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2827        else:           #OK  ellipsoidal crystallites
2828            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2829            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2830            H = np.array(refl[:3])
2831            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2832            Sgam = const/lenR
2833            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2834                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2835                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2836        gamDict[phfx+'Size;mx'] = Sgam  #OK
2837        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
2838               
2839        #microstrain derivatives               
2840        if calcControls[phfx+'MustrainType'] == 'isotropic':
2841            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2842            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
2843            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
2844        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2845            H = np.array(refl[:3])
2846            P = np.array(calcControls[phfx+'MustrainAxis'])
2847            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2848            Si = parmDict[phfx+'Mustrain;i']
2849            Sa = parmDict[phfx+'Mustrain;a']
2850            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
2851            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2852            Mgam = gami/sqtrm
2853            dsi = -gami*Si*cosP**2/sqtrm**3
2854            dsa = -gami*Sa*sinP**2/sqtrm**3
2855            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2856            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2857            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2858            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2859        else:       #generalized - P.W. Stephens model OK
2860            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2861            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
2862            Sum = 0
2863            for i,strm in enumerate(Strms):
2864                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2865                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2866                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2867            Mgam = const*np.sqrt(Sum)
2868            for i in range(len(Strms)):
2869                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2870                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2       
2871        gamDict[phfx+'Mustrain;mx'] = Mgam
2872        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2873       
2874    return sigDict,gamDict
2875       
2876def GetReflPos(refl,im,wave,A,pfx,hfx,phfx,calcControls,parmDict):
2877    'Needs a doc string'
2878    if im:
2879        h,k,l,m = refl[:4]
2880        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2881        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
2882    else:
2883        h,k,l = refl[:3]
2884        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
2885    refl[4+im] = d
2886    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:
2887        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2888        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2889        if 'Bragg' in calcControls[hfx+'instType']:
2890            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2891                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2892        else:               #Debye-Scherrer - simple but maybe not right
2893            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+(parmDict[hfx+'DisplaceY']+parmDict[phfx+'LayerDisp'])*sind(pos))
2894    elif 'T' in calcControls[hfx+'histType']:
2895        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
2896        #do I need sample position effects - maybe?
2897    return pos
2898
2899def GetReflPosDerv(refl,im,wave,A,pfx,hfx,phfx,calcControls,parmDict):
2900    'Needs a doc string'
2901    dpr = 180./np.pi
2902    if im:
2903        h,k,l,m = refl[:4]
2904        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2905        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
2906        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
2907    else:
2908        m = 0
2909        h,k,l = refl[:3]       
2910        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2911    dst = np.sqrt(dstsq)
2912    dsp = 1./dst
2913    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:
2914        pos = refl[5+im]-parmDict[hfx+'Zero']
2915        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2916        dpdw = const*dst
2917        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
2918        dpdZ = 1.0
2919        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2920            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
2921        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2922        if 'Bragg' in calcControls[hfx+'instType']:
2923            dpdSh = -4.*shft*cosd(pos/2.0)
2924            dpdTr = -shft*sind(pos)*100.0
2925            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
2926        else:               #Debye-Scherrer - simple but maybe not right
2927            dpdXd = -shft*cosd(pos)
2928            dpdYd = -shft*sind(pos)
2929            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
2930    elif 'T' in calcControls[hfx+'histType']:
2931        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
2932        dpdZ = 1.0
2933        dpdDC = dsp
2934        dpdDA = dsp**2
2935        dpdDB = 1./dsp
2936        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2937            2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2.
2938        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2939           
2940def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2941    'Needs a doc string'
2942    laue = SGData['SGLaue']
2943    uniq = SGData['SGUniq']
2944    h,k,l = refl[:3]
2945    if laue in ['m3','m3m']:
2946        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2947            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2948    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2949        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2950    elif laue in ['3R','3mR']:
2951        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2952    elif laue in ['4/m','4/mmm']:
2953        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2954    elif laue in ['mmm']:
2955        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2956    elif laue in ['2/m']:
2957        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2958        if uniq == 'a':
2959            Dij += parmDict[phfx+'D23']*k*l
2960        elif uniq == 'b':
2961            Dij += parmDict[phfx+'D13']*h*l
2962        elif uniq == 'c':
2963            Dij += parmDict[phfx+'D12']*h*k
2964    else:
2965        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2966            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2967    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:
2968        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2969    else:
2970        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
2971           
2972def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2973    'Needs a doc string'
2974    laue = SGData['SGLaue']
2975    uniq = SGData['SGUniq']
2976    h,k,l = refl[:3]
2977    if laue in ['m3','m3m']:
2978        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2979            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2980    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2981        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2982    elif laue in ['3R','3mR']:
2983        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2984    elif laue in ['4/m','4/mmm']:
2985        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2986    elif laue in ['mmm']:
2987        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2988    elif laue in ['2/m']:
2989        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2990        if uniq == 'a':
2991            dDijDict[phfx+'D23'] = k*l
2992        elif uniq == 'b':
2993            dDijDict[phfx+'D13'] = h*l
2994        elif uniq == 'c':
2995            dDijDict[phfx+'D12'] = h*k
2996    else:
2997        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2998            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2999    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:
3000        for item in dDijDict:
3001            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
3002    else:
3003        for item in dDijDict:
3004            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
3005    return dDijDict
3006   
3007def GetDij(phfx,SGData,parmDict):
3008    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
3009    return G2spc.HStrainVals(HSvals,SGData)
3010               
3011def GetFobsSq(Histograms,Phases,parmDict,calcControls):
3012    '''Compute the observed structure factors for Powder histograms and store in reflection array
3013    Multiprocessing support added
3014    '''
3015    if GSASIIpath.GetConfigValue('Show_timing',False):
3016        starttime = time.time() #; print 'start GetFobsSq'
3017    histoList = list(Histograms.keys())
3018    histoList.sort()
3019    Ka2 = shl = lamRatio = kRatio = None
3020    for histogram in histoList:
3021        if 'PWDR' in histogram[:4]:
3022            Histogram = Histograms[histogram]
3023            hId = Histogram['hId']
3024            hfx = ':%d:'%(hId)
3025            Limits = calcControls[hfx+'Limits']
3026            if 'C' in calcControls[hfx+'histType']:
3027                shl = max(parmDict[hfx+'SH/L'],0.0005)
3028                Ka2 = False
3029                kRatio = 0.0
3030                if hfx+'Lam1' in list(parmDict.keys()):
3031                    Ka2 = True
3032                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3033                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
3034            x,y,w,yc,yb,yd = Histogram['Data']
3035            xMask = ma.getmaskarray(x)
3036            xB = np.searchsorted(x,Limits[0])
3037            xF = np.searchsorted(x,Limits[1])
3038            ymb = np.array(y-yb)
3039            ymb = np.where(ymb,ymb,1.0)
3040            ycmb = np.array(yc-yb)
3041            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
3042            refLists = Histogram['Reflection Lists']
3043            for phase in refLists:
3044                if phase not in Phases:     #skips deleted or renamed phases silently!
3045                    continue
3046                Phase = Phases[phase]
3047                if histogram not in Phase['Histograms']:
3048                    continue
3049                im = 0
3050                if Phase['General'].get('Modulated',False):
3051                    im = 1
3052                pId = Phase['pId']
3053                phfx = '%d:%d:'%(pId,hId)
3054                refDict = refLists[phase]
3055                sumFo = 0.0
3056                sumdF = 0.0
3057                sumFosq = 0.0
3058                sumdFsq = 0.0
3059                sumInt = 0.0
3060                nExcl = 0
3061                # test to see if we are using multiprocessing below
3062                useMP,ncores = G2mp.InitMP()
3063                if len(refDict['RefList']) < 100: useMP = False       
3064                if useMP: # multiprocessing: create a set of initialized Python processes
3065                    MPpool = mp.Pool(G2mp.ncores,G2mp.InitFobsSqGlobals,
3066                                    [x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2])
3067                    profArgs = [[] for i in range(G2mp.ncores)]
3068                else:
3069                    G2mp.InitFobsSqGlobals(x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2)
3070                if 'C' in calcControls[hfx+'histType']:
3071                    # are we multiprocessing?
3072                    for iref,refl in enumerate(refDict['RefList']):
3073                        if useMP: 
3074                            profArgs[iref%G2mp.ncores].append((refl,iref))
3075                        else:
3076                            icod = G2mp.ComputeFobsSqCW(refl,iref)
3077                            if type(icod) is tuple:
3078                                refl[8+im] = icod[0]
3079                                sumInt += icod[1]
3080                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3081                            elif icod == -1:
3082                                refl[3+im] *= -1
3083                                nExcl += 1
3084                            elif icod == -2:
3085                                break
3086                    if useMP:
3087                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):
3088                            sumInt += sInt
3089                            for refl8im,irefl in resList:
3090                                if refl8im is None:
3091                                    refDict['RefList'][irefl][3+im] *= -1
3092                                    nExcl += 1
3093                                else:
3094                                    refDict['RefList'][irefl][8+im] = refl8im
3095                                    if parmDict[phfx+'LeBail']:
3096                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3097                elif 'T' in calcControls[hfx+'histType']:
3098                    for iref,refl in enumerate(refDict['RefList']):
3099                        if useMP: 
3100                            profArgs[iref%G2mp.ncores].append((refl,iref))
3101                        else:
3102                            icod = G2mp.ComputeFobsSqTOF(refl,iref)
3103                            if type(icod) is tuple:
3104                                refl[8+im] = icod[0]
3105                                sumInt += icod[1]
3106                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3107                            elif icod == -1:
3108                                refl[3+im] *= -1
3109                                nExcl += 1
3110                            elif icod == -2:
3111                                break
3112                    if useMP:
3113                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):
3114                            sumInt += sInt
3115                            for refl8im,irefl in resList:
3116                                if refl8im is None:
3117                                    refDict['RefList'][irefl][3+im] *= -1
3118                                    nExcl += 1
3119                                else:
3120                                    refDict['RefList'][irefl][8+im] = refl8im
3121                                    if parmDict[phfx+'LeBail']:
3122                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3123                elif 'B' in calcControls[hfx+'histType']:
3124                    for iref,refl in enumerate(refDict['RefList']):
3125                        if useMP: 
3126                            profArgs[iref%G2mp.ncores].append((refl,iref))
3127                        else:
3128                            icod = G2mp.ComputeFobsSqPink(refl,iref)
3129                            if type(icod) is tuple:
3130                                refl[8+im] = icod[0]
3131                                sumInt += icod[1]
3132                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3133                            elif icod == -1:
3134                                refl[3+im] *= -1
3135                                nExcl += 1
3136                            elif icod == -2:
3137                                break
3138                    if useMP:
3139                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqPinkbatch,profArgs):
3140                            sumInt += sInt
3141                            for refl8im,irefl in resList:
3142                                if refl8im is None:
3143                                    refDict['RefList'][irefl][3+im] *= -1
3144                                    nExcl += 1
3145                                else:
3146                                    refDict['RefList'][irefl][8+im] = refl8im
3147                                    if parmDict[phfx+'LeBail']:
3148                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3149                if useMP: MPpool.terminate()
3150                sumFo = 0.0
3151                sumdF = 0.0
3152                sumFosq = 0.0
3153                sumdFsq = 0.0
3154                for iref,refl in enumerate(refDict['RefList']):
3155                    Fo = np.sqrt(np.abs(refl[8+im]))
3156                    Fc = np.sqrt(np.abs(refl[9+im]))
3157                    sumFo += Fo
3158                    sumFosq += refl[8+im]**2
3159                    sumdF += np.abs(Fo-Fc)
3160                    sumdFsq += (refl[8+im]-refl[9+im])**2
3161                if sumFo:
3162                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3163                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3164                else:
3165                    Histogram['Residuals'][phfx+'Rf'] = 100.
3166                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
3167                Histogram['Residuals'][phfx+'sumInt'] = sumInt
3168                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])-nExcl
3169                Histogram['Residuals']['hId'] = hId
3170        elif 'HKLF' in histogram[:4]:
3171            Histogram = Histograms[histogram]
3172            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
3173    if GSASIIpath.GetConfigValue('Show_timing',False):
3174        print ('GetFobsSq t=',time.time()-starttime)
3175               
3176def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup,histogram=None):
3177    'Computes the powder pattern for a histogram based on contributions from all used phases'
3178    if GSASIIpath.GetConfigValue('Show_timing',False): starttime = time.time()
3179   
3180    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
3181        U = parmDict[hfx+'U']
3182        V = parmDict[hfx+'V']
3183        W = parmDict[hfx+'W']
3184        X = parmDict[hfx+'X']
3185        Y = parmDict[hfx+'Y']
3186        Z = parmDict[hfx+'Z']
3187        tanPos = tand(refl[5+im]/2.0)
3188        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3189        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3190        sig = max(0.001,sig)
3191        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam+Z     #save peak gamma
3192        gam = max(0.001,gam)
3193        return sig,gam
3194               
3195    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
3196        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
3197            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']*refl[4+im]
3198        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2+parmDict[hfx+'Z']
3199        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3200        sig += Ssig
3201        gam += Sgam
3202        return sig,gam
3203       
3204    def GetReflAlpBet(refl,im,hfx,parmDict):
3205        alp = parmDict[hfx+'alpha']/refl[4+im]
3206        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
3207        return alp,bet
3208       
3209    def GetPinkReflAlpBet(refl,im,hfx,parmDict):
3210        tanPos = tand(refl[5+im]/2.0)
3211        alp = max(0.1,parmDict[hfx+'alpha-0']+parmDict[hfx+'alpha-1']*tanPos)
3212        bet = max(0.001,parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']*tanPos)
3213        return alp,bet
3214
3215    hId = Histogram['hId']
3216    hfx = ':%d:'%(hId)
3217    bakType = calcControls[hfx+'bakType']
3218    fixback = Histogram['Background'][1].get('fixback',None)
3219    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x,fixback)
3220    yc = np.zeros_like(yb)
3221    # cw = np.diff(ma.getdata(x))
3222    # cw = np.append(cw,cw[-1])
3223       
3224    if 'C' in calcControls[hfx+'histType']:   
3225        shl = max(parmDict[hfx+'SH/L'],0.002)
3226        Ka2 = False
3227        if hfx+'Lam1' in (parmDict.keys()):
3228            wave = parmDict[hfx+'Lam1']
3229            Ka2 = True
3230            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3231            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3232        else:
3233            wave = parmDict[hfx+'Lam']
3234    elif 'B' in calcControls[hfx+'histType']:
3235            wave = parmDict[hfx+'Lam']
3236    else:
3237        shl = 0.
3238    for phase in Histogram['Reflection Lists']:
3239        refDict = Histogram['Reflection Lists'][phase]
3240        if phase not in Phases:     #skips deleted or renamed phases silently!
3241            continue
3242        Phase = Phases[phase]
3243        if histogram and not histogram in Phase['Histograms']:
3244            continue
3245        pId = Phase['pId']
3246        pfx = '%d::'%(pId)
3247        phfx = '%d:%d:'%(pId,hId)
3248        hfx = ':%d:'%(hId)
3249        SGData = Phase['General']['SGData']
3250        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3251        im = 0
3252        if Phase['General'].get('Modulated',False):
3253            SSGData = Phase['General']['SSGData']
3254            im = 1  #offset in SS reflection list
3255            #??
3256        Dij = GetDij(phfx,SGData,parmDict)
3257        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]  #TODO: need to do something if Dij << 0.
3258        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3259        if np.any(np.diag(G)<0.):
3260            msg = 'Invalid metric tensor for phase #{}\n   ({})'.format(
3261                pId,Phase['General']['Name'])
3262        elif np.any(np.isnan(A)):
3263            msg = 'Infinite metric tensor for phase #{}\n   ({})'.format(
3264                pId,Phase['General']['Name'])
3265        else:
3266            msg = None
3267        if msg:
3268            print('\nInvalid cell metric tensor for phase #{} ({})\n'.format(
3269                pId,Phase['General']['Name']),
3270                'values (A): {:.3f} {:.3f} {:.3f} {:.3f} {:.3f} {:.3f}\n'
3271                .format(*A))
3272            raise G2obj.G2Exception('Error: '+msg+
3273                        ' See console.\nCheck for refinement of conflicting variables')
3274        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matrices
3275        Vst = np.sqrt(nl.det(G))    #V*
3276        if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
3277            if im:
3278                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
3279            elif parmDict[pfx+'isMag'] and 'N' in calcControls[hfx+'histType']:
3280                MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3281            else:
3282                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3283        badPeak = False
3284        # test to see if we are using multiprocessing here
3285        useMP,ncores = G2mp.InitMP()
3286        if len(refDict['RefList']) < 100: useMP = False       
3287        if useMP: # multiprocessing: create a set of initialized Python processes
3288            MPpool = mp.Pool(ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
3289            profArgs = [[] for i in range(ncores)]
3290        if 'C' in calcControls[hfx+'histType']:
3291            for iref,refl in enumerate(refDict['RefList']):
3292                if im:
3293                    h,k,l,m = refl[:4]
3294                else:
3295                    h,k,l = refl[:3]
3296                Uniq = np.inner(refl[:3],SGMT)
3297                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,phfx,calcControls,parmDict)         #corrected reflection position
3298                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
3299                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3300                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3301                refl[11+im] *= Vst*Lorenz
3302                 
3303                if Phase['General'].get('doPawley'):
3304                    try:
3305                        if im:
3306                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3307                        else:
3308                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3309                        refl[9+im] = parmDict[pInd]
3310                    except KeyError:
3311#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3312                        continue
3313                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
3314                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3315                iFin = np.searchsorted(x,refl[5+im]+fmax)
3316                if not iBeg+iFin:       #peak below low limit - skip peak
3317                    continue
3318                elif not iBeg-iFin:     #peak above high limit - done
3319                    break
3320                elif iBeg > iFin:   #bad peak coeff - skip
3321                    badPeak = True
3322                    continue
3323                if useMP:
3324                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin,1.))
3325                else:
3326                    fp = G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))[0]
3327                    yc[iBeg:iFin] += refl[11+im]*refl[9+im]*fp   #>90% of time spent here
3328                if Ka2:
3329                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3330                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
3331                    iBeg = np.searchsorted(x,pos2-fmin)
3332                    iFin = np.searchsorted(x,pos2+fmax)
3333                    if not iBeg+iFin:       #peak below low limit - skip peak
3334                        continue
3335                    elif not iBeg-iFin:     #peak above high limit - done
3336                        return yc,yb
3337                    elif iBeg > iFin: