source: trunk/GSASIIstrMath.py @ 4869

Last change on this file since 4869 was 4869, checked in by toby, 7 months ago

.Raise() is causing crashes, use only on wx.Dialog or wx.Frame

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