source: trunk/GSASIIstrMath.py @ 5274

Last change on this file since 5274 was 5274, checked in by toby, 19 months ago

implement save of intensities by phase (partials)

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