source: trunk/GSASIIstrMath.py @ 4840

Last change on this file since 4840 was 4840, checked in by toby, 9 months ago

LeBail? issues: bug on zero cycles fixed; new LeBail? fit menu item; HessianLSQ: remove extra func eval on 0 cycles; fix noprint options; misc cleanups

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