source: trunk/GSASIIstrMath.py @ 4531

Last change on this file since 4531 was 4531, checked in by vondreele, 17 months ago

remove q-plot & d-plot options for Limits PWDR plot
some cleanup in mag inncom. sf cal

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 225.4 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2020-07-21 17:41:16 +0000 (Tue, 21 Jul 2020) $
8# $Author: vondreele $
9# $Revision: 4531 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 4531 2020-07-21 17:41:16Z vondreele $
12########### SVN repository information ###################
13from __future__ import division, print_function
14import time
15import copy
16import numpy as np
17import numpy.ma as ma
18import numpy.linalg as nl
19import scipy.stats as st
20import multiprocessing as mp
21import GSASIIpath
22GSASIIpath.SetVersionNumber("$Revision: 4531 $")
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 (xyz & 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            TLS = RBObj['ThermalMotion']
95            if 'T' in TLS[0]:
96                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
97                    TLS[1][i] = parmDict[pfx+pt+rbsx]
98            if 'L' in TLS[0]:
99                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
100                    TLS[1][i+6] = parmDict[pfx+pt+rbsx]
101            if 'S' in TLS[0]:
102                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
103                    TLS[1][i+12] = parmDict[pfx+pt+rbsx]
104            if 'U' in TLS[0]:
105                TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
106            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
107            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
108            for i,x in enumerate(XYZ):
109                atId = RBObj['Ids'][i]
110                for j in [0,1,2]:
111                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
112                if UIJ[i][0] == 'A':
113                    for j in range(6):
114                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
115                elif UIJ[i][0] == 'I':
116                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
117           
118        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
119            jrb = RRBIds.index(RBObj['RBId'])
120            rbsx = str(irb)+':'+str(jrb)
121            for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
122                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
123            for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
124                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]               
125            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
126            TLS = RBObj['ThermalMotion']
127            if 'T' in TLS[0]:
128                for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
129                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
130            if 'L' in TLS[0]:
131                for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
132                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
133            if 'S' in TLS[0]:
134                for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
135                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
136            if 'U' in TLS[0]:
137                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
138            for itors,tors in enumerate(RBObj['Torsions']):
139                tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
140            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
141            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
142            for i,x in enumerate(XYZ):
143                atId = RBObj['Ids'][i]
144                for j in [0,1,2]:
145                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
146                if UIJ[i][0] == 'A':
147                    for j in range(6):
148                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
149                elif UIJ[i][0] == 'I':
150                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
151                   
152def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
153    'Needs a doc string'
154    atxIds = ['dAx:','dAy:','dAz:']
155    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
156    OIds = ['Oa:','Oi:','Oj:','Ok:']
157    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
158    if not RBIds['Vector'] and not RBIds['Residue']:
159        return
160    VRBIds = RBIds['Vector']
161    RRBIds = RBIds['Residue']
162    RBData = rigidbodyDict
163    for item in parmDict:
164        if 'RB' in item:
165            dFdvDict[item] = 0.        #NB: this is a vector which is no. refl. long & must be filled!
166    General = Phase['General']
167    cx,ct,cs,cia = General['AtomPtrs']
168    cell = General['Cell'][1:7]
169    Amat,Bmat = G2lat.cell2AB(cell)
170    rpd = np.pi/180.
171    rpd2 = rpd**2
172    g = nl.inv(np.inner(Bmat,Bmat))
173    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
174        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
175    AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
176    pfx = str(Phase['pId'])+'::'
177    RBModels =  Phase['RBModels']
178    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
179        VModel = RBData['Vector'][RBObj['RBId']]
180        Q = RBObj['Orient'][0]
181        jrb = VRBIds.index(RBObj['RBId'])
182        rbsx = str(irb)+':'+str(jrb)
183        dXdv = []
184        for iv in range(len(VModel['VectMag'])):
185            dCdv = []
186            for vec in VModel['rbVect'][iv]:
187                dCdv.append(G2mth.prodQVQ(Q,vec))
188            dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
189        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
190        for ia,atId in enumerate(RBObj['Ids']):
191            atNum = AtLookup[atId]
192            dx = 0.00001
193            for iv in range(len(VModel['VectMag'])):
194                for ix in [0,1,2]:
195                    dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
196            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
197                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
198            for iv in range(4):
199                Q[iv] -= dx
200                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
201                Q[iv] += 2.*dx
202                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
203                Q[iv] -= dx
204                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
205                for ix in [0,1,2]:
206                    dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
207            X = G2mth.prodQVQ(Q,Cart[ia])
208            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
209            dFdu = G2lat.U6toUij(dFdu.T)
210            dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))           
211            dFdu = G2lat.UijtoU6(dFdu)
212            atNum = AtLookup[atId]
213            if 'T' in RBObj['ThermalMotion'][0]:
214                for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
215                    dFdvDict[pfx+name+rbsx] += dFdu[i]
216            if 'L' in RBObj['ThermalMotion'][0]:
217                dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
218                dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
219                dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
220                dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
221                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
222                dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2*(-dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
223                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
224                dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2*(-dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
225                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
226            if 'S' in RBObj['ThermalMotion'][0]:
227                dFdvDict[pfx+'RBVS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
228                dFdvDict[pfx+'RBVS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
229                dFdvDict[pfx+'RBVS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
230                dFdvDict[pfx+'RBVS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
231                dFdvDict[pfx+'RBVS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
232                dFdvDict[pfx+'RBVS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
233                dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
234                dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
235            if 'U' in RBObj['ThermalMotion'][0]:
236                dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
237
238
239    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
240        Q = RBObj['Orient'][0]
241        jrb = RRBIds.index(RBObj['RBId'])
242        torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
243        rbsx = str(irb)+':'+str(jrb)
244        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
245        for itors,tors in enumerate(RBObj['Torsions']):     #derivative error?
246            tname = pfx+'RBRTr;'+str(itors)+':'+rbsx           
247            orId,pvId = torData[itors][:2]
248            pivotVec = Cart[orId]-Cart[pvId]
249            QA = G2mth.AVdeg2Q(-0.001,pivotVec)
250            QB = G2mth.AVdeg2Q(0.001,pivotVec)
251            for ir in torData[itors][3]:
252                atNum = AtLookup[RBObj['Ids'][ir]]
253                rVec = Cart[ir]-Cart[pvId]
254                dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
255                dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
256                for ix in [0,1,2]:
257                    dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
258        for ia,atId in enumerate(RBObj['Ids']):
259            atNum = AtLookup[atId]
260            dx = 0.00001
261            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
262                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
263            for iv in range(4):
264                Q[iv] -= dx
265                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
266                Q[iv] += 2.*dx
267                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
268                Q[iv] -= dx
269                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
270                for ix in [0,1,2]:
271                    dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
272            X = G2mth.prodQVQ(Q,Cart[ia])
273            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
274            dFdu = G2lat.U6toUij(dFdu.T)
275            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
276            dFdu = G2lat.UijtoU6(dFdu)
277            atNum = AtLookup[atId]
278            if 'T' in RBObj['ThermalMotion'][0]:
279                for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
280                    dFdvDict[pfx+name+rbsx] += dFdu[i]
281            if 'L' in RBObj['ThermalMotion'][0]:
282                dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
283                dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
284                dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
285                dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
286                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
287                dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2*(dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
288                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
289                dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2*(dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
290                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
291            if 'S' in RBObj['ThermalMotion'][0]:
292                dFdvDict[pfx+'RBRS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
293                dFdvDict[pfx+'RBRS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
294                dFdvDict[pfx+'RBRS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
295                dFdvDict[pfx+'RBRS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
296                dFdvDict[pfx+'RBRS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
297                dFdvDict[pfx+'RBRS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
298                dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
299                dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
300            if 'U' in RBObj['ThermalMotion'][0]:
301                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
302   
303################################################################################
304##### Penalty & restraint functions
305################################################################################
306
307def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
308    'Compute user-supplied and built-in restraint functions'
309    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
310    pNames = []
311    pVals = []
312    pWt = []
313    negWt = {}
314    pWsum = {}
315    pWnum = {}
316    for phase in Phases:
317        pId = Phases[phase]['pId']
318        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
319        General = Phases[phase]['General']
320        cx,ct,cs,cia = General['AtomPtrs']
321        textureData = General['SH Texture']
322        SGData = General['SGData']
323        Atoms = Phases[phase]['Atoms']
324        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
325        cell = General['Cell'][1:7]
326        Amat,Bmat = G2lat.cell2AB(cell)
327        if phase not in restraintDict:
328            continue
329        phaseRest = restraintDict[phase]
330        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
331            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
332            ['ChemComp','Sites'],['Texture','HKLs'],['General','General'],]
333        for name,rest in names:
334            pWsum[name] = 0.
335            pWnum[name] = 0
336            if name not in phaseRest:
337                continue
338            itemRest = phaseRest[name]
339            if itemRest[rest] and itemRest['Use']:
340                wt = itemRest.get('wtFactor',1.)
341                if name in ['Bond','Angle','Plane','Chiral']:
342                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
343                        pNames.append(str(pId)+':'+name+':'+str(i))
344                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
345                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
346                        if name == 'Bond':
347                            calc = G2mth.getRestDist(XYZ,Amat)
348                        elif name == 'Angle':
349                            calc = G2mth.getRestAngle(XYZ,Amat)
350                        elif name == 'Plane':
351                            calc = G2mth.getRestPlane(XYZ,Amat)
352                        elif name == 'Chiral':
353                            calc = G2mth.getRestChiral(XYZ,Amat)
354                        pVals.append(obs-calc)
355                        pWt.append(wt/esd**2)
356                        pWsum[name] += wt*((obs-calc)/esd)**2
357                        pWnum[name] += 1
358                elif name in ['Torsion','Rama']:
359                    coeffDict = itemRest['Coeff']
360                    for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
361                        pNames.append(str(pId)+':'+name+':'+str(i))
362                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
363                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
364                        if name == 'Torsion':
365                            tor = G2mth.getRestTorsion(XYZ,Amat)
366                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
367                        else:
368                            phi,psi = G2mth.getRestRama(XYZ,Amat)
369                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
370                        pVals.append(restr)
371                        pWt.append(wt/esd**2)
372                        pWsum[name] += wt*(restr/esd)**2
373                        pWnum[name] += 1
374                elif name == 'ChemComp':
375                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
376                        pNames.append(str(pId)+':'+name+':'+str(i))
377                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
378                        frac = np.array(G2mth.GetAtomFracByID(pId,parmDict,AtLookup,indx))
379                        calc = np.sum(mul*frac*factors)
380                        pVals.append(obs-calc)
381                        pWt.append(wt/esd**2)                   
382                        pWsum[name] += wt*((obs-calc)/esd)**2
383                        pWnum[name] += 1
384                elif name == 'Texture':
385                    SHkeys = list(textureData['SH Coeff'][1].keys())
386                    SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
387                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
388                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
389                    for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
390                        PH = np.array(hkl)
391                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
392                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
393                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
394                        Z1 = ma.masked_greater(Z,0.0)           #is this + or -?
395                        IndZ1 = np.array(ma.nonzero(Z1))
396                        for ind in IndZ1.T:
397                            pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
398                            pVals.append(Z1[ind[0]][ind[1]])
399                            pWt.append(wt/esd1**2)
400                            pWsum[name] += wt*(-Z1[ind[0]][ind[1]]/esd1)**2
401                            pWnum[name] += 1
402                        if ifesd2:
403                            Z2 = 1.-Z
404                            for ind in np.ndindex(grid,grid):
405                                pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
406                                pVals.append(Z2[ind[0]][ind[1]])
407                                pWt.append(wt/esd2**2)
408                                pWsum[name] += wt*(Z2/esd2)**2
409                                pWnum[name] += 1
410                elif name == 'General':
411                    for i,(eq,obs,esd) in enumerate(itemRest[rest]):
412                        pNames.append(str(pId)+':'+name+':'+str(i))
413                        calcobj = G2obj.ExpressionCalcObj(eq)
414                        calcobj.SetupCalc(parmDict)
415                        calc = calcobj.EvalExpression()
416                        pVals.append(obs-calc)
417                        pWt.append(wt/esd**2)                   
418                        pWsum[name] += wt*((obs-calc)/esd)**2
419                        pWnum[name] += 1
420       
421    for phase in Phases:
422        name = 'SH-Pref.Ori.'
423        pId = Phases[phase]['pId']
424        General = Phases[phase]['General']
425        SGData = General['SGData']
426        cell = General['Cell'][1:7]
427        pWsum[name] = 0.0
428        pWnum[name] = 0
429        for hist in Phases[phase]['Histograms']:
430            if not Phases[phase]['Histograms'][hist]['Use']:
431                continue
432            if hist in Histograms and 'PWDR' in hist:
433                hId = Histograms[hist]['hId']
434                phfx = '%d:%d:'%(pId,hId)
435                if calcControls[phfx+'poType'] == 'SH':
436                    toler = calcControls[phfx+'SHtoler']
437                    wt = 1./toler**2
438                    HKLs = np.array(calcControls[phfx+'SHhkl'])
439                    SHnames = calcControls[phfx+'SHnames']
440                    SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
441                    for i,PH in enumerate(HKLs):
442                        phi,beta = G2lat.CrsAng(PH,cell,SGData)
443                        SH3Coef = {}
444                        for item in SHcof:
445                            L,N = eval(item.strip('C'))
446                            SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
447                        ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
448                        X = np.linspace(0,90.0,26)
449                        Y = ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0)       #+ or -?
450                        IndY = ma.nonzero(Y)
451                        for ind in IndY[0]:
452                            pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind]))
453                            pVals.append(Y[ind])
454                            pWt.append(wt)
455                            pWsum[name] += wt*(Y[ind])**2
456                            pWnum[name] += 1
457    pWsum['PWLref'] = 0.
458    pWnum['PWLref'] = 0
459    for item in varyList:
460        if 'PWLref' in item and parmDict[item] < 0.:
461            pId = int(item.split(':')[0])
462            if negWt[pId]:
463                pNames.append(item)
464                pVals.append(parmDict[item])
465                pWt.append(negWt[pId])
466                pWsum['PWLref'] += negWt[pId]*(parmDict[item])**2
467                pWnum['PWLref'] += 1
468    pVals = np.array(pVals)
469    pWt = np.array(pWt)         #should this be np.sqrt?
470    return pNames,pVals,pWt,pWsum,pWnum
471   
472def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
473    '''Compute derivatives on user-supplied and built-in restraint
474    (penalty) functions
475
476    where pNames is list of restraint labels
477
478    returns pDerv with partial derivatives by variable# in varList and
479       restraint# in pNames (pDerv[variable#][restraint#])
480    '''
481    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
482    pDerv = np.zeros((len(varyList),len(pVal)))
483    for pName in pNames: # loop over restraints
484        if 'General' == pName.split(':')[1]:
485            # initialize for General restraint(s) here
486            GeneralInit = True
487            parmDict0 = parmDict.copy()
488            # setup steps for each parameter
489            stepDict = {}
490            for parm in varyList:
491                stepDict[parm] = G2obj.getVarStep(parm,parmDict)
492            break
493    for phase in Phases:
494#        if phase not in restraintDict:
495#            continue
496        pId = Phases[phase]['pId']
497        General = Phases[phase]['General']
498        cx,ct,cs,cia = General['AtomPtrs']
499        SGData = General['SGData']
500        Atoms = Phases[phase]['Atoms']
501        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
502        cell = General['Cell'][1:7]
503        Amat,Bmat = G2lat.cell2AB(cell)
504        textureData = General['SH Texture']
505
506        SHkeys = list(textureData['SH Coeff'][1].keys())
507        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
508        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
509        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
510        sam = SamSym[textureData['Model']]
511        phaseRest = restraintDict.get(phase,{})
512        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
513            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
514            'ChemComp':'Sites','Texture':'HKLs'}
515        lasthkl = np.array([0,0,0])
516        for ip,pName in enumerate(pNames): # loop over restraints
517            pnames = pName.split(':')
518            if pId == int(pnames[0]):
519                name = pnames[1]
520                if 'PWL' in pName:
521                    pDerv[varyList.index(pName)][ip] += 1.
522                    continue
523                elif 'SH-' in pName:
524                    continue
525                Id = int(pnames[2]) 
526                itemRest = phaseRest[name]
527                if name in ['Bond','Angle','Plane','Chiral']:
528                    indx,ops,obs,esd = itemRest[names[name]][Id]
529                    dNames = []
530                    for ind in indx:
531                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
532                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
533                    if name == 'Bond':
534                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
535                    elif name == 'Angle':
536                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
537                    elif name == 'Plane':
538                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
539                    elif name == 'Chiral':
540                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
541                elif name in ['Torsion','Rama']:
542                    coffDict = itemRest['Coeff']
543                    indx,ops,cofName,esd = itemRest[names[name]][Id]
544                    dNames = []
545                    for ind in indx:
546                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
547                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
548                    if name == 'Torsion':
549                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
550                    else:
551                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
552                elif name == 'ChemComp':
553                    indx,factors,obs,esd = itemRest[names[name]][Id]
554                    dNames = []
555                    for ind in indx:
556                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
557                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
558                        deriv = mul*factors
559                elif 'Texture' in name:
560                    deriv = []
561                    dNames = []
562                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][Id]
563                    hkl = np.array(hkl)
564                    if np.any(lasthkl-hkl):
565                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
566                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
567                        lasthkl = copy.copy(hkl)                       
568                    if 'unit' in name:
569                        pass
570                    else:
571                        gam = float(pnames[3])
572                        psi = float(pnames[4])
573                        for SHname in ODFln:
574                            l,m,n = eval(SHname[1:])
575                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
576                            dNames += [str(pId)+'::'+SHname]
577                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
578                elif name == 'General':
579                    deriv = []
580                    dNames = []
581                    eq,obs,esd = itemRest[name][Id]
582                    calcobj = G2obj.ExpressionCalcObj(eq)
583                    parmlist = list(eq.assgnVars.values()) # parameters used in this expression
584                    for parm in parmlist: # expand list if any parms are determined by constraints
585                        if parm in G2mv.dependentVars:
586                            parmlist += G2mv.independentVars
587                            break
588                    for ind,var in enumerate(varyList):
589                        drv = 0
590                        if var in parmlist:
591                            step = stepDict.get(var,1e-5)
592                            calc = []
593                            # apply step to parameter
594                            oneparm = True
595                            for s in -step,2*step:
596                                parmDict[var] += s
597                                # extend shift if needed to other parameters
598                                if var in G2mv.independentVars:
599                                    G2mv.Dict2Map(parmDict,[])
600                                    oneparm = False
601                                elif var in G2mv.dependentVars:
602                                    G2mv.Map2Dict(parmDict,[])
603                                    oneparm = False
604                                if 'RB' in var:
605                                    ApplyRBModels(parmDict,Phases,rigidbodyDict)
606# test
607                                    oneparm = False
608                                calcobj.SetupCalc(parmDict)
609                                calc.append(calcobj.EvalExpression())
610                            drv = (calc[1]-calc[0])*.5/step
611                            # restore the dict
612                            if oneparm:
613                                parmDict[var] = parmDict0[var]
614                            else:
615                                parmDict = parmDict0.copy()
616                        else:
617                            drv = 0
618                        pDerv[ind][ip] = drv
619                # Add derivatives into matrix, if needed
620                for dName,drv in zip(dNames,deriv):
621                    try:
622                        ind = varyList.index(dName)
623                        pDerv[ind][ip] += drv
624                    except ValueError:
625                        pass
626       
627        lasthkl = np.array([0,0,0])
628        for ip,pName in enumerate(pNames):
629            deriv = []
630            dNames = []
631            pnames = pName.split(':')
632            if 'SH-' in pName and pId == int(pnames[0]):
633                hId = int(pnames[1])
634                phfx = '%d:%d:'%(pId,hId)
635                psi = float(pnames[4])
636                HKLs = calcControls[phfx+'SHhkl']
637                SHnames = calcControls[phfx+'SHnames']
638                SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
639                hkl = np.array(HKLs[int(pnames[3])])     
640                if np.any(lasthkl-hkl):
641                    phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
642                    SH3Coef = {}
643                    for item in SHcof:
644                        L,N = eval(item.strip('C'))
645                        SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
646                    ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
647                    lasthkl = copy.copy(hkl)                       
648                for SHname in SHnames:
649                    l,n = eval(SHname[1:])
650                    SH3name = 'C%d,0,%d'%(l,n)
651                    Ksl = G2lat.GetKsl(l,0,'0',psi,0.0)[0]
652                    dNames += [phfx+SHname]
653                    deriv.append(ODFln[SH3name][0]*Ksl/SHcof[SHname])
654            for dName,drv in zip(dNames,deriv):
655                try:
656                    ind = varyList.index(dName)
657                    pDerv[ind][ip] += drv
658                except ValueError:
659                    pass
660    return pDerv
661
662################################################################################
663##### Function & derivative calculations
664################################################################################       
665                   
666def GetAtomFXU(pfx,calcControls,parmDict):
667    'Needs a doc string'
668    Natoms = calcControls['Natoms'][pfx]
669    Tdata = Natoms*[' ',]
670    Mdata = np.zeros(Natoms)
671    IAdata = Natoms*[' ',]
672    Fdata = np.zeros(Natoms)
673    Xdata = np.zeros((3,Natoms))
674    dXdata = np.zeros((3,Natoms))
675    Uisodata = np.zeros(Natoms)
676    Uijdata = np.zeros((6,Natoms))
677    Gdata = np.zeros((3,Natoms))
678    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
679        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
680        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
681        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
682        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5],
683        'AMx:':Gdata[0],'AMy:':Gdata[1],'AMz:':Gdata[2],}
684    for iatm in range(Natoms):
685        for key in keys:
686            parm = pfx+key+str(iatm)
687            if parm in parmDict:
688                keys[key][iatm] = parmDict[parm]
689    Fdata = np.where(Fdata,Fdata,1.e-8)         #avoid divide by zero in derivative calc.
690    Gdata = np.where(Gdata,Gdata,1.e-8)         #avoid divide by zero in derivative calc.
691   
692    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata
693   
694def GetAtomSSFXU(pfx,calcControls,parmDict):
695    'Needs a doc string'
696    Natoms = calcControls['Natoms'][pfx]
697    maxSSwave = calcControls['maxSSwave'][pfx]
698    Nwave = {'F':maxSSwave['Sfrac'],'X':maxSSwave['Spos'],'Y':maxSSwave['Spos'],'Z':maxSSwave['Spos'],
699        'U':maxSSwave['Sadp'],'M':maxSSwave['Smag'],'T':maxSSwave['Spos']}
700    XSSdata = np.zeros((6,maxSSwave['Spos'],Natoms))
701    FSSdata = np.zeros((2,maxSSwave['Sfrac'],Natoms))
702    USSdata = np.zeros((12,maxSSwave['Sadp'],Natoms))
703    MSSdata = np.zeros((6,maxSSwave['Smag'],Natoms))
704    waveTypes = []
705    keys = {'Fsin:':FSSdata[0],'Fcos:':FSSdata[1],'Fzero:':FSSdata[0],'Fwid:':FSSdata[1],
706        'Tmin:':XSSdata[0],'Tmax:':XSSdata[1],'Xmax:':XSSdata[2],'Ymax:':XSSdata[3],'Zmax:':XSSdata[4],
707        'Xsin:':XSSdata[0],'Ysin:':XSSdata[1],'Zsin:':XSSdata[2],'Xcos:':XSSdata[3],'Ycos:':XSSdata[4],'Zcos:':XSSdata[5],
708        'U11sin:':USSdata[0],'U22sin:':USSdata[1],'U33sin:':USSdata[2],'U12sin:':USSdata[3],'U13sin:':USSdata[4],'U23sin:':USSdata[5],
709        'U11cos:':USSdata[6],'U22cos:':USSdata[7],'U33cos:':USSdata[8],'U12cos:':USSdata[9],'U13cos:':USSdata[10],'U23cos:':USSdata[11],
710        'MXsin:':MSSdata[0],'MYsin:':MSSdata[1],'MZsin:':MSSdata[2],'MXcos:':MSSdata[3],'MYcos:':MSSdata[4],'MZcos:':MSSdata[5]}
711    for iatm in range(Natoms):
712        wavetype = [parmDict.get(pfx+kind+'waveType:'+str(iatm),'') for kind in ['F','P','A','M']]
713        waveTypes.append(wavetype)
714        for key in keys:
715            for m in range(Nwave[key[0]]):
716                parm = pfx+key+str(iatm)+':%d'%(m)
717                if parm in parmDict:
718                    keys[key][m][iatm] = parmDict[parm]
719    return waveTypes,FSSdata,XSSdata,USSdata,MSSdata
720   
721def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
722    ''' Compute structure factors for all h,k,l for phase
723    puts the result, F^2, in each ref[8] in refList
724    operates on blocks of 100 reflections for speed
725    input:
726   
727    :param dict refDict: where
728        'RefList' list where each ref = h,k,l,it,d,...
729        'FF' dict of form factors - filed in below
730    :param np.array G:      reciprocal metric tensor
731    :param str pfx:    phase id string
732    :param dict SGData: space group info. dictionary output from SpcGroup
733    :param dict calcControls:
734    :param dict ParmDict:
735
736    '''       
737    phfx = pfx.split(':')[0]+hfx
738    ast = np.sqrt(np.diag(G))
739    Mast = twopisq*np.multiply.outer(ast,ast)
740    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
741    SGT = np.array([ops[1] for ops in SGData['SGOps']])
742    FFtables = calcControls['FFtables']
743    BLtables = calcControls['BLtables']
744    Amat,Bmat = G2lat.Gmat2AB(G)
745    Flack = 1.0
746    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
747        Flack = 1.-2.*parmDict[phfx+'Flack']
748    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
749    TwDict = refDict.get('TwDict',{})           
750    if 'S' in calcControls[hfx+'histType']:
751        NTL = calcControls[phfx+'NTL']
752        NM = calcControls[phfx+'TwinNMN']+1
753        TwinLaw = calcControls[phfx+'TwinLaw']
754        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
755        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
756    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
757        GetAtomFXU(pfx,calcControls,parmDict)
758    if not Xdata.size:          #no atoms in phase!
759        return
760    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
761        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
762    elif 'X' in calcControls[hfx+'histType']:
763        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
764        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
765    Uij = np.array(G2lat.U6toUij(Uijdata))
766    bij = Mast*Uij.T
767    blkSize = 100       #no. of reflections in a block - size seems optimal
768    nRef = refDict['RefList'].shape[0]
769    SQ = 1./(2.*refDict['RefList'].T[4])**2
770    if 'N' in calcControls[hfx+'histType']:
771        dat = G2el.getBLvalues(BLtables)
772        refDict['FF']['El'] = list(dat.keys())
773        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
774    else:       #'X'
775        dat = G2el.getFFvalues(FFtables,0.)
776        refDict['FF']['El'] = list(dat.keys())
777        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
778        for iel,El in enumerate(refDict['FF']['El']):
779            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
780#reflection processing begins here - big arrays!
781    iBeg = 0
782    while iBeg < nRef:
783        iFin = min(iBeg+blkSize,nRef)
784        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
785        H = refl.T[:3]                          #array(blkSize,3)
786        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
787        TwMask = np.any(H,axis=-1)
788        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
789            for ir in range(blkSize):
790                iref = ir+iBeg
791                if iref in TwDict:
792                    for i in TwDict[iref]:
793                        for n in range(NTL):
794                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
795            TwMask = np.any(H,axis=-1)
796        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
797        SQfactor = 4.0*SQ*twopisq               #ditto prev.
798        if 'T' in calcControls[hfx+'histType']:
799            if 'P' in calcControls[hfx+'histType']:
800                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
801            else:
802                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
803            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
804            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
805        Uniq = np.inner(H,SGMT)
806        Phi = np.inner(H,SGT)
807        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
808        sinp = np.sin(phase)
809        cosp = np.cos(phase)
810        biso = -SQfactor*Uisodata[:,nxs]
811        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
812        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
813        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
814        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
815        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
816        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
817        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
818        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
819            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
820            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
821        else:
822            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
823            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
824        fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
825        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
826        if SGData['SGInv']: #centrosymmetric; B=0
827            fbs[0] *= 0.
828            fas[1] *= 0.
829        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
830            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
831            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
832        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
833            if len(TwinLaw) > 1:
834                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
835                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
836                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
837                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
838            else:   # checked correct!!
839                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
840                refl.T[7] = np.copy(refl.T[9])               
841                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
842#                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
843        iBeg += blkSize
844#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
845   
846def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
847    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
848    faster than StructureFactorDerv - correct for powders/nontwins!!
849    input:
850   
851    :param dict refDict: where
852        'RefList' list where each ref = h,k,l,it,d,...
853        'FF' dict of form factors - filled in below
854    :param np.array G:      reciprocal metric tensor
855    :param str hfx:    histogram id string
856    :param str pfx:    phase id string
857    :param dict SGData: space group info. dictionary output from SpcGroup
858    :param dict calcControls:
859    :param dict parmDict:
860   
861    :returns: dict dFdvDict: dictionary of derivatives
862    '''
863    phfx = pfx.split(':')[0]+hfx
864    ast = np.sqrt(np.diag(G))
865    Mast = twopisq*np.multiply.outer(ast,ast)
866    SGMT = np.array([ops[0] for ops in SGData['SGOps']])
867    SGT = np.array([ops[1] for ops in SGData['SGOps']])
868    FFtables = calcControls['FFtables']
869    BLtables = calcControls['BLtables']
870    Amat,Bmat = G2lat.Gmat2AB(G)
871    nRef = len(refDict['RefList'])
872    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
873        GetAtomFXU(pfx,calcControls,parmDict)
874    if not Xdata.size:          #no atoms in phase!
875        return {}
876    mSize = len(Mdata)
877    FF = np.zeros(len(Tdata))
878    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
879        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
880    elif 'X' in calcControls[hfx+'histType']:
881        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
882        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
883    Uij = np.array(G2lat.U6toUij(Uijdata))
884    bij = Mast*Uij.T
885    dFdvDict = {}
886    dFdfr = np.zeros((nRef,mSize))
887    dFdx = np.zeros((nRef,mSize,3))
888    dFdui = np.zeros((nRef,mSize))
889    dFdua = np.zeros((nRef,mSize,6))
890    dFdbab = np.zeros((nRef,2))
891    dFdfl = np.zeros((nRef))
892    Flack = 1.0
893    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
894        Flack = 1.-2.*parmDict[phfx+'Flack']
895#    time0 = time.time()
896#reflection processing begins here - big arrays!
897    iBeg = 0
898    blkSize = 32       #no. of reflections in a block - optimized for speed
899    while iBeg < nRef:
900        iFin = min(iBeg+blkSize,nRef)
901        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
902        H = refl.T[:3].T
903        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
904        SQfactor = 8.0*SQ*np.pi**2
905        if 'T' in calcControls[hfx+'histType']:
906            if 'P' in calcControls[hfx+'histType']:
907                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
908            else:
909                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
910            FP = np.repeat(FP.T,len(SGT),axis=0)
911            FPP = np.repeat(FPP.T,len(SGT),axis=0)
912        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
913        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
914        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
915        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
916        Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
917        Phi = np.inner(H,SGT)
918        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
919        sinp = np.sin(phase)        #refBlk x nOps x nAtoms
920        cosp = np.cos(phase)
921        occ = Mdata*Fdata/len(SGT)
922        biso = -SQfactor*Uisodata[:,nxs]
923        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
924        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
925        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
926        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
927        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])      #Nref*Nops,3,3
928        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))     #Nref,Nops,6
929        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
930        if len(FPP.shape) > 1:
931            fotp = np.reshape(FPP,cosp.shape)*Tcorr
932        else:
933            fotp = FPP*Tcorr     
934        if 'T' in calcControls[hfx+'histType']:
935            fa = np.array([fot*cosp,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
936            fb = np.array([fot*sinp,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
937        else:
938            fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr])
939            fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr])
940        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,refBlk,nTwins)
941        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
942        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,refBlk,nEqv,nAtoms)
943        fbx = np.array([fot*cosp,-fotp*sinp])
944        #sum below is over Uniq
945        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,nAtom) Fdata != 0 avoids /0. problem
946        dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
947        dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
948        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nAtoms)
949        dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
950        # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
951        if not SGData['SGInv']:
952            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
953            dfbdba = np.sum(-sinp*Tcorr,axis=-2)
954            dfadfl = np.sum(np.sum(-fotp*sinp,axis=-1),axis=-1)
955            dfbdfl = np.sum(np.sum(fotp*cosp,axis=-1),axis=-1)
956            dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
957            dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
958            dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2)
959        else:
960            dfbdfr = np.zeros_like(dfadfr)
961            dfbdx = np.zeros_like(dfadx)
962            dfbdui = np.zeros_like(dfadui)
963            dfbdua = np.zeros_like(dfadua)
964            dfbdba = np.zeros_like(dfadba)
965            dfadfl = 0.0
966            dfbdfl = 0.0
967        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
968        SA = fas[0]+fas[1]
969        SB = fbs[0]+fbs[1]
970        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
971            dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/len(SGMT)
972            dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0)
973            dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0)
974            dFdua[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadua+fbs[:,:,nxs,nxs]*dfbdua,axis=0)
975        else:
976            dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
977            dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
978            dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
979            dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])
980            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
981        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)])+ \
982                            fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
983        iBeg += blkSize
984#    print 'derv time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
985        #loop over atoms - each dict entry is list of derivatives for all the reflections
986    for i in range(len(Mdata)):
987        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
988        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
989        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
990        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
991        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
992        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
993        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
994        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
995        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
996        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
997        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
998    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
999    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1000    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1001    return dFdvDict
1002   
1003def MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1004    ''' Compute neutron magnetic structure factors for all h,k,l for phase
1005    puts the result, F^2, in each ref[8] in refList
1006    operates on blocks of 100 reflections for speed
1007    input:
1008   
1009    :param dict refDict: where
1010        'RefList' list where each ref = h,k,l,it,d,...
1011        'FF' dict of form factors - filed in below
1012    :param np.array G:      reciprocal metric tensor
1013    :param str pfx:    phase id string
1014    :param dict SGData: space group info. dictionary output from SpcGroup
1015    :param dict calcControls:
1016    :param dict ParmDict:
1017       
1018    :returns: copy of new refList - used in calculating numerical derivatives
1019
1020    '''       
1021    g = nl.inv(G)
1022    ast = np.sqrt(np.diag(G))
1023    ainv = np.sqrt(np.diag(g))
1024    GS = G/np.outer(ast,ast)
1025    Ginv = g/np.outer(ainv,ainv)
1026    uAmat = G2lat.Gmat2AB(GS)[0]
1027    Mast = twopisq*np.multiply.outer(ast,ast)
1028    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1029    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1030    Ncen = len(SGData['SGCen'])
1031    Nops = len(SGMT)*Ncen
1032    if not SGData['SGFixed']:
1033        Nops *= (1+SGData['SGInv'])
1034    MFtables = calcControls['MFtables']
1035    Bmat = G2lat.Gmat2AB(G)[1]
1036    TwinLaw = np.ones(1)
1037#    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
1038#    TwDict = refDict.get('TwDict',{})           
1039#    if 'S' in calcControls[hfx+'histType']:
1040#        NTL = calcControls[phfx+'NTL']
1041#        NM = calcControls[phfx+'TwinNMN']+1
1042#        TwinLaw = calcControls[phfx+'TwinLaw']
1043#        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1044#        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1045    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1046        GetAtomFXU(pfx,calcControls,parmDict)
1047    if not Xdata.size:          #no atoms in phase!
1048        return
1049    Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T])
1050    Gdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2)).T            #apply sym. ops.
1051    if SGData['SGInv'] and not SGData['SGFixed']:
1052        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
1053    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering--> [Mxyz,nops,natms]
1054    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
1055    Mag = np.tile(Mag[:,nxs],Nops).#make Mag same length as Gdata
1056    Kdata = np.inner(Gdata.T,uAmat).T
1057    Kmean = np.mean(np.sqrt(np.sum(Kdata**2,axis=0)),axis=0)
1058    Kdata /= Kmean     #Cartesian unit vectors
1059    Uij = np.array(G2lat.U6toUij(Uijdata))
1060    bij = Mast*Uij.T
1061    blkSize = 100       #no. of reflections in a block - size seems optimal
1062    nRef = refDict['RefList'].shape[0]
1063    SQ = 1./(2.*refDict['RefList'].T[4])**2
1064    refDict['FF']['El'] = list(MFtables.keys())
1065    refDict['FF']['MF'] = np.zeros((nRef,len(MFtables)))
1066    for iel,El in enumerate(refDict['FF']['El']):
1067        refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1068#reflection processing begins here - big arrays!
1069    iBeg = 0
1070    while iBeg < nRef:
1071        iFin = min(iBeg+blkSize,nRef)
1072        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1073        H = refl.T[:3].T                          #array(blkSize,3)
1074#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
1075#        TwMask = np.any(H,axis=-1)
1076#        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1077#            for ir in range(blkSize):
1078#                iref = ir+iBeg
1079#                if iref in TwDict:
1080#                    for i in TwDict[iref]:
1081#                        for n in range(NTL):
1082#                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1083#            TwMask = np.any(H,axis=-1)
1084        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
1085        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1086        Uniq = np.inner(H,SGMT)
1087        Phi = np.inner(H,SGT)
1088        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1089        biso = -SQfactor*Uisodata[:,nxs]
1090        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
1091        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1092        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
1093        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1094        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1095        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
1096        if SGData['SGInv']:
1097            if not SGData['SGFixed']:
1098                mphase = np.hstack((phase,-phase))  #OK
1099            else:
1100                mphase = phase
1101        else:
1102            mphase = phase                    #
1103        mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
1104        mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
1105        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
1106        cosm = np.cos(mphase)                               #ditto
1107        HM = np.inner(Bmat,H)                             #put into cartesian space
1108        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Kdata = MAGS & HM = UVEC in magstrfc.for both OK
1109        eDotK = np.sum(HM[:,:,nxs,nxs]*Kdata[:,nxs,:,:],axis=0)
1110        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
1111        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
1112        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
1113        fams = np.sum(np.sum(fam,axis=-1),axis=-1)     #Mxyz,Nref  Sum(sum(fam,atoms),ops)
1114        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)     #ditto
1115        refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)   #Sum(fams**2,Mxyz) Re + Im
1116        refl.T[7] = np.copy(refl.T[9])               
1117        refl.T[10] = atan2d(fbms[0],fams[0]) #- what is phase for mag refl?
1118#        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
1119#            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
1120#            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1121#        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
1122#            if len(TwinLaw) > 1:
1123#                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
1124#                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
1125#                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
1126#                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
1127#            else:   # checked correct!!
1128#                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
1129#                refl.T[7] = np.copy(refl.T[9])               
1130#                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1131##                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
1132        iBeg += blkSize
1133#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
1134    return copy.deepcopy(refDict['RefList'])
1135
1136def MagStructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1137    '''Compute magnetic structure factor derivatives numerically - for powders/nontwins only
1138    input:
1139   
1140    :param dict refDict: where
1141        'RefList' list where each ref = h,k,l,it,d,...
1142        'FF' dict of form factors - filled in below
1143    :param np.array G:      reciprocal metric tensor
1144    :param str hfx:    histogram id string
1145    :param str pfx:    phase id string
1146    :param dict SGData: space group info. dictionary output from SpcGroup
1147    :param dict calcControls:
1148    :param dict parmDict:
1149   
1150    :returns: dict dFdvDict: dictionary of magnetic derivatives
1151    '''
1152   
1153    trefDict = copy.deepcopy(refDict)
1154    dM = 1.e-6
1155    dFdvDict = {}
1156    for parm in parmDict:
1157        if 'AM' in parm:
1158            parmDict[parm] += dM
1159            prefList = MagStructureFactor2(trefDict,G,hfx,pfx,SGData,calcControls,parmDict)
1160            parmDict[parm] -= 2*dM
1161            mrefList = MagStructureFactor2(trefDict,G,hfx,pfx,SGData,calcControls,parmDict)
1162            parmDict[parm] += dM
1163            dFdvDict[parm] = (prefList[:,9]-mrefList[:,9])/(2.*dM)
1164    return dFdvDict
1165           
1166def MagStructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1167    '''Compute nonmagnetic structure factor derivatives on blocks of reflections in magnetic structures - for powders/nontwins only
1168    input:
1169   
1170    :param dict refDict: where
1171        'RefList' list where each ref = h,k,l,it,d,...
1172        'FF' dict of form factors - filled in below
1173    :param np.array G:      reciprocal metric tensor
1174    :param str hfx:    histogram id string
1175    :param str pfx:    phase id string
1176    :param dict SGData: space group info. dictionary output from SpcGroup
1177    :param dict calcControls:
1178    :param dict parmDict:
1179   
1180    :returns: dict dFdvDict: dictionary of derivatives
1181    '''
1182   
1183    g = nl.inv(G)
1184    ast = np.sqrt(np.diag(G))
1185    ainv = np.sqrt(np.diag(g))
1186    GS = G/np.outer(ast,ast)
1187    Ginv = g/np.outer(ainv,ainv)
1188    uAmat = G2lat.Gmat2AB(GS)[0]
1189    Mast = twopisq*np.multiply.outer(ast,ast)
1190    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1191    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1192    Ncen = len(SGData['SGCen'])
1193    Nops = len(SGMT)*Ncen
1194    if not SGData['SGFixed']:
1195        Nops *= (1+SGData['SGInv'])
1196    Bmat = G2lat.Gmat2AB(G)[1]
1197    nRef = len(refDict['RefList'])
1198    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1199        GetAtomFXU(pfx,calcControls,parmDict)
1200    if not Xdata.size:          #no atoms in phase!
1201        return {}
1202    mSize = len(Mdata)
1203    Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T])
1204    Gones = np.ones_like(Gdata)
1205    Gdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2)).T            #apply sym. ops.
1206    Gones = np.inner(Gones.T,SGMT).T
1207    if SGData['SGInv'] and not SGData['SGFixed']:
1208        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
1209        Gones = np.hstack((Gones,-Gones))       #inversion if any
1210    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
1211    Gones = np.hstack([Gones for icen in range(Ncen)])        #dup over cell centering
1212    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip
1213    Gones = SGData['MagMom'][nxs,:,nxs]*Gones   #flip vectors according to spin flip
1214    Mag = np.tile(Mag[:,nxs],Nops).#make Mag same length as Gdata
1215    Kdata = np.inner(Gdata.T,uAmat).T     #Cartesian unit vectors
1216    Kmean = np.mean(np.sqrt(np.sum(Kdata**2,axis=0)),axis=0)
1217    Kdata /= Kmean
1218    Uij = np.array(G2lat.U6toUij(Uijdata))
1219    bij = Mast*Uij.T
1220    dFdvDict = {}
1221    dFdfr = np.zeros((nRef,mSize))
1222    dFdx = np.zeros((nRef,mSize,3))
1223    dFdui = np.zeros((nRef,mSize))
1224    dFdua = np.zeros((nRef,mSize,6))
1225    time0 = time.time()
1226#reflection processing begins here - big arrays!
1227    iBeg = 0
1228    blkSize = 5       #no. of reflections in a block - optimized for speed
1229    while iBeg < nRef:
1230        iFin = min(iBeg+blkSize,nRef)
1231        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1232        H = refl.T[:3].T
1233        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1234        SQfactor = 8.0*SQ*np.pi**2
1235        Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
1236        Phi = np.inner(H,SGT)
1237        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1238        occ = Mdata*Fdata/Nops
1239        biso = -SQfactor*Uisodata[:,nxs]
1240        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
1241        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1242        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
1243        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1244        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))
1245        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1246        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1247        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
1248        if SGData['SGInv']:
1249            if not SGData['SGFixed']:
1250                mphase = np.hstack((phase,-phase))  #OK
1251                Uniq = np.hstack((Uniq,-Uniq))      #Nref,Nops,hkl
1252                Hij = np.hstack((Hij,Hij))
1253            else:
1254                mphase = phase
1255        else:
1256            mphase = phase                    #
1257        Hij = np.concatenate(np.array([Hij for cen in SGData['SGCen']]),axis=1)
1258        Uniq = np.hstack([Uniq for cen in SGData['SGCen']])
1259        mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
1260        mphase = np.concatenate(mphase,axis=1)              #Nref,Nop,Natm
1261        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
1262        cosm = np.cos(mphase)                               #ditto
1263        HM = np.inner(Bmat.T,H)                             #put into cartesian space
1264        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #unit cartesian vector for H
1265        eDotK = np.sum(HM[:,:,nxs,nxs]*Kdata[:,nxs,:,:],axis=0)
1266        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK
1267       
1268        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #Mxyz,Nref,Nop,Natm
1269        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]
1270        fams = np.sum(np.sum(fam,axis=-1),axis=-1)                      #Mxyz,Nref
1271        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)
1272        famx = -Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*sinm[nxs,:,:,:]   #Mxyz,Nref,Nops,Natom
1273        fbmx = Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*cosm[nxs,:,:,:]
1274        #sums below are over Nops - real part
1275        dfadfr = np.sum(fam/occ,axis=2)        #array(Mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem deriv OK
1276        dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*famx[:,:,:,:,nxs],axis=2)          #deriv OK
1277        dfadui = np.sum(-SQfactor[:,nxs,nxs]*fam,axis=2) #array(Ops,refBlk,nAtoms)  deriv OK
1278        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)            #deriv OK
1279        # imaginary part; array(3,refBlk,nAtom,3) & array(3,refBlk,nAtom,6)
1280        dfbdfr = np.sum(fbm/occ,axis=2)        #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem
1281        dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*fbmx[:,:,:,:,nxs],axis=2)
1282        dfbdui = np.sum(-SQfactor[:,nxs,nxs]*fbm,axis=2) #array(Ops,refBlk,nAtoms)
1283        dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbm[:,:,:,:,nxs],axis=2)
1284        #accumulate derivatives   
1285        dFdfr[iBeg:iFin] = 2.*np.sum((fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/Nops,axis=0) #ok
1286        dFdx[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadx+fbms[:,:,nxs,nxs]*dfbdx,axis=0)         #ok
1287        dFdui[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui,axis=0)              #ok
1288        dFdua[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua,axis=0)      #ok
1289        iBeg += blkSize
1290    print (' %d derivative time %.4f\r'%(nRef,time.time()-time0))
1291        #loop over atoms - each dict entry is list of derivatives for all the reflections
1292    for i in range(len(Mdata)):
1293        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1294        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1295        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1296        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1297        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1298        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1299        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1300        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1301        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1302        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1303        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1304    return dFdvDict
1305       
1306def StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1307    '''Compute structure factor derivatives on blocks of reflections - for twins only
1308    faster than StructureFactorDervTw
1309    input:
1310   
1311    :param dict refDict: where
1312        'RefList' list where each ref = h,k,l,it,d,...
1313        'FF' dict of form factors - filled in below
1314    :param np.array G:      reciprocal metric tensor
1315    :param str hfx:    histogram id string
1316    :param str pfx:    phase id string
1317    :param dict SGData: space group info. dictionary output from SpcGroup
1318    :param dict calcControls:
1319    :param dict parmDict:
1320   
1321    :returns: dict dFdvDict: dictionary of derivatives
1322    '''
1323    phfx = pfx.split(':')[0]+hfx
1324    ast = np.sqrt(np.diag(G))
1325    Mast = twopisq*np.multiply.outer(ast,ast)
1326    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1327    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1328    FFtables = calcControls['FFtables']
1329    BLtables = calcControls['BLtables']
1330    TwDict = refDict.get('TwDict',{})           
1331    NTL = calcControls[phfx+'NTL']
1332    NM = calcControls[phfx+'TwinNMN']+1
1333    TwinLaw = calcControls[phfx+'TwinLaw']
1334    TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1335    TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1336    nTwin = len(TwinLaw)       
1337    nRef = len(refDict['RefList'])
1338    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1339        GetAtomFXU(pfx,calcControls,parmDict)
1340    if not Xdata.size:          #no atoms in phase!
1341        return {}
1342    mSize = len(Mdata)
1343    FF = np.zeros(len(Tdata))
1344    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
1345        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1346    elif 'X' in calcControls[hfx+'histType']:
1347        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1348        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1349    Uij = np.array(G2lat.U6toUij(Uijdata))
1350    bij = Mast*Uij.T
1351    dFdvDict = {}
1352    dFdfr = np.zeros((nRef,nTwin,mSize))
1353    dFdx = np.zeros((nRef,nTwin,mSize,3))
1354    dFdui = np.zeros((nRef,nTwin,mSize))
1355    dFdua = np.zeros((nRef,nTwin,mSize,6))
1356    dFdbab = np.zeros((nRef,nTwin,2))
1357    dFdtw = np.zeros((nRef,nTwin))
1358    time0 = time.time()
1359#reflection processing begins here - big arrays!
1360    iBeg = 0
1361    blkSize = 16       #no. of reflections in a block - optimized for speed
1362    while iBeg < nRef:
1363        iFin = min(iBeg+blkSize,nRef)
1364        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1365        H = refl.T[:3]
1366        H = np.inner(H.T,TwinLaw)   #array(3,nTwins)
1367        TwMask = np.any(H,axis=-1)
1368        for ir in range(blkSize):
1369            iref = ir+iBeg
1370            if iref in TwDict:
1371                for i in TwDict[iref]:
1372                    for n in range(NTL):
1373                        H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1374        TwMask = np.any(H,axis=-1)
1375        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1376        SQfactor = 8.0*SQ*np.pi**2
1377        if 'T' in calcControls[hfx+'histType']:
1378            if 'P' in calcControls[hfx+'histType']:
1379                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1380            else:
1381                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1382            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
1383            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
1384        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1385        Bab = np.repeat(parmDict[phfx+'BabA']*dBabdA,len(SGT)*nTwin)
1386        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1387        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
1388        Uniq = np.inner(H,SGMT)             # (nTwin,nSGOp,3)
1389        Phi = np.inner(H,SGT)
1390        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1391        sinp = np.sin(phase)
1392        cosp = np.cos(phase)
1393        occ = Mdata*Fdata/len(SGT)
1394        biso = -SQfactor*Uisodata[:,nxs]
1395        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
1396        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1397        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1398        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,nTwin,len(SGT),6))
1399        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1400        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*Mdata*Fdata/len(SGMT)
1401        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
1402        fotp = FPP*Tcorr       
1403        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
1404            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(FPP,sinp.shape)*sinp*Tcorr])
1405            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(FPP,cosp.shape)*cosp*Tcorr])
1406        else:
1407            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-FPP*sinp*Tcorr])
1408            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,FPP*cosp*Tcorr])
1409        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
1410        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
1411        if SGData['SGInv']: #centrosymmetric; B=0
1412            fbs[0] *= 0.
1413            fas[1] *= 0.
1414        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,nRef,ntwi,nEqv,nAtoms)
1415        fbx = np.array([fot*cosp,-fotp*sinp])
1416        #sum below is over Uniq
1417        dfadfr = np.sum(np.sum(fa/occ,axis=-2),axis=0)        #array(2,nRef,ntwin,nAtom) Fdata != 0 avoids /0. problem
1418        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1419        dfadui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),axis=0)           
1420        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'
1421        dfadua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fa[:,:,:,:,:,nxs],axis=-3),axis=0) 
1422        if not SGData['SGInv']:
1423            dfbdfr = np.sum(np.sum(fb/occ,axis=-2),axis=0)        #Fdata != 0 avoids /0. problem
1424            dfadba /= 2.
1425#            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
1426            dfbdui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),axis=0)
1427            dfbdx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fbx[:,:,:,:,:,nxs],axis=-3),axis=0)
1428            dfbdua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fb[:,:,:,:,:,nxs],axis=-3),axis=0)
1429        else:
1430            dfbdfr = np.zeros_like(dfadfr)
1431            dfbdx = np.zeros_like(dfadx)
1432            dfbdui = np.zeros_like(dfadui)
1433            dfbdua = np.zeros_like(dfadua)
1434#            dfbdba = np.zeros_like(dfadba)
1435        SA = fas[0]+fas[1]
1436        SB = fbs[0]+fbs[1]
1437        dFdfr[iBeg:iFin] = ((2.*TwMask*SA)[:,:,nxs]*dfadfr+(2.*TwMask*SB)[:,:,nxs]*dfbdfr)*Mdata[nxs,nxs,:]/len(SGMT)
1438        dFdx[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadx+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdx
1439        dFdui[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs]*dfadui+(2.*TwMask*SB)[:,:,nxs]*dfbdui
1440        dFdua[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadua+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdua
1441        if SGData['SGInv']: #centrosymmetric; B=0
1442            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
1443        else:               
1444            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
1445#        dFdbab[iBeg:iFin] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1446#            fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1447        iBeg += blkSize
1448    print (' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0))
1449    #loop over atoms - each dict entry is list of derivatives for all the reflections
1450    for i in range(len(Mdata)):     #these all OK
1451        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
1452        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
1453        dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
1454        dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
1455        dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
1456        dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
1457        dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
1458        dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
1459        dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
1460        dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
1461        dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
1462    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1463    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1464    for i in range(nTwin):
1465        dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
1466    return dFdvDict
1467   
1468def SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1469    '''
1470    Compute super structure factors for all h,k,l,m for phase - no twins
1471    puts the result, F^2, in each ref[9] in refList
1472    works on blocks of 32 reflections for speed
1473    input:
1474   
1475    :param dict refDict: where
1476        'RefList' list where each ref = h,k,l,m,it,d,...
1477        'FF' dict of form factors - filed in below
1478    :param np.array G:      reciprocal metric tensor
1479    :param str pfx:    phase id string
1480    :param dict SGData: space group info. dictionary output from SpcGroup
1481    :param dict calcControls:
1482    :param dict ParmDict:
1483
1484    '''
1485    phfx = pfx.split(':')[0]+hfx
1486    ast = np.sqrt(np.diag(G))
1487    GS = G/np.outer(ast,ast)
1488    uAmat,uBmat = G2lat.Gmat2AB(GS)
1489    Mast = twopisq*np.multiply.outer(ast,ast)   
1490    SGInv = SGData['SGInv']
1491    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1492    Nops = len(SGMT)*(1+SGData['SGInv'])
1493    if SGData['SGGray']:
1494        Nops *= 2
1495    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1496    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1497    SSCen = SSGData['SSGCen']
1498    FFtables = calcControls['FFtables']
1499    BLtables = calcControls['BLtables']
1500    MFtables = calcControls['MFtables']
1501    Amat,Bmat = G2lat.Gmat2AB(G)
1502    Flack = 1.0
1503    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1504        Flack = 1.-2.*parmDict[phfx+'Flack']
1505    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1506        GetAtomFXU(pfx,calcControls,parmDict)
1507    if not Xdata.size:          #no atoms in phase!
1508        return
1509    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1510    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)      #NB: Mmod is ReIm,Mxyz,Ntau,Natm
1511    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1512
1513    if parmDict[pfx+'isMag']:       #This part correct for making modulated mag moments on equiv atoms - Mmod matched drawing & Bilbao drawings
1514   
1515        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
1516        if SGData['SGGray']:
1517            mXYZ = np.hstack((mXYZ,mXYZ))
1518# this done in wrong place - should be
1519        MmodA,MmodB = G2mth.MagMod(glTau,mXYZ,modQ,MSSdata,SGData,SSGData)  #Ntau,Nops,Natm,Mxyz cos,sim parts sum matches drawing
1520        Mmod = MmodA+MmodB
1521       
1522        if not SGData['SGGray']:    #for fixed Mx,My,Mz
1523            GSdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2))  #apply sym. ops.--> Natm,Nops,Nxyz
1524            if SGData['SGInv'] and not SGData['SGFixed']:   #inversion if any
1525                GSdata = np.hstack((GSdata,-GSdata))     
1526            GSdata = np.hstack([GSdata for cen in SSCen])        #dup over cell centering - Natm,Nops,Mxyz
1527            GSdata = SGData['MagMom'][nxs,:,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
1528            GSdata = np.swapaxes(GSdata,0,1)    #Nop,Natm,Mxyz
1529            Mmod += GSdata[nxs,:,:,:]
1530
1531#        Kdata = np.inner(Mmod,uAmat)          #Ntau,Nops,Natm,Mxyz
1532#        Kmean = np.mean(np.sqrt(np.sum(Kdata**2,axis=0)),axis=0)
1533#        Kdata /= Kmean     #Ntau,Nops,Natm,Mxyz  Cartesian unit vectors along mag moments
1534       
1535    FF = np.zeros(len(Tdata))
1536    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
1537        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1538    elif 'X' in calcControls[hfx+'histType']:
1539        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1540        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1541    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1542    bij = Mast*Uij
1543    blkSize = 48       #no. of reflections in a block
1544    nRef = refDict['RefList'].shape[0]
1545    SQ = 1./(2.*refDict['RefList'].T[5])**2
1546    if 'N' in calcControls[hfx+'histType']:
1547        dat = G2el.getBLvalues(BLtables)
1548        refDict['FF']['El'] = list(dat.keys())
1549        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1550        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1551        for iel,El in enumerate(refDict['FF']['El']):
1552            if El in MFtables:
1553                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1554    else:
1555        dat = G2el.getFFvalues(FFtables,0.)
1556        refDict['FF']['El'] = list(dat.keys())
1557        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1558        for iel,El in enumerate(refDict['FF']['El']):
1559            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1560#    time0 = time.time()
1561#reflection processing begins here - big arrays!
1562    iBeg = 0
1563    while iBeg < nRef:
1564        iFin = min(iBeg+blkSize,nRef)
1565        mRef= iFin-iBeg
1566        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1567        H = refl.T[:4]                          #array(blkSize,4)
1568        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1569        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1570        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1571        Uniq = np.inner(H.T,SSGMT)
1572        UniqP = np.inner(HP.T,SGMT)
1573        Phi = np.inner(H.T,SSGT)
1574        if SGInv and not SGData['SGFixed']:   #if centro - expand HKL sets
1575            Uniq = np.hstack((Uniq,-Uniq))
1576            Phi = np.hstack((Phi,-Phi))
1577            UniqP = np.hstack((UniqP,-UniqP))
1578        if 'T' in calcControls[hfx+'histType']:
1579            if 'P' in calcControls[hfx+'histType']:
1580                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1581            else:
1582                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1583            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
1584            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
1585        Bab = 0.
1586        if phfx+'BabA' in parmDict:
1587            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
1588        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1589        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
1590        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1591        phase = np.hstack([phase for cen in SSCen])
1592        sinp = np.sin(phase)
1593        cosp = np.cos(phase)
1594        biso = -SQfactor*Uisodata[:,nxs]
1595        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T
1596        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1597        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1598        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1599
1600        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:           
1601           
1602            phasem = twopi*np.inner(mXYZ,HP.T).T    #2pi(Q.r)
1603            cosm = np.cos(phasem)                   #Nref,nops,natm
1604            sinm = np.sin(phasem)
1605            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1606            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Mdata*Fdata*MF/(2*Nops)     #Nref,Natm
1607            HM = np.inner(uAmat,HP.T)                    #put into cartesian space X||H,Z||H*L; uAmat better than uBmat
1608            eM = (HM/np.sqrt(np.sum(HM**2,axis=0))).T    # normalize  HP  Nref,hkl=Unit vectors || Q
1609
1610            fam0 = 0.
1611            fbm0 = 0.
1612            if not SGData['SGGray']:     #correct -fixed Mx,My,Mz contribution             
1613                fam0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*cosm[:,:,:,nxs]    #Nref,Nops,Natm,Mxyz
1614                fbm0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*sinm[:,:,:,nxs]
1615#  calc mag. structure factors; Nref,Ntau,Nops,Natm,Mxyz                           
1616            fams = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*cosm[i,nxs,:,:,nxs]+   
1617                H[3,i]*MmodB*sinm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
1618                       
1619            fbms = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*sinm[i,nxs,:,:,nxs]+   
1620                H[3,i]*MmodB*cosm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
1621
1622            if not SGData['SGGray']:           
1623                fams += fam0[:,nxs,:,:,:]
1624                fbms += fbm0[:,nxs,:,:,:]
1625               
1626# this is the best, but not exactly right
1627#sum ops & atms                               
1628            fasm = np.sum(np.sum(fams,axis=-2),axis=-2)    #Nref,Ntau,Mxyz; sum ops & atoms
1629            fbsm = np.sum(np.sum(fbms,axis=-2),axis=-2)
1630#put into cartesian space
1631            facm = np.inner(fasm,uBmat.T)
1632            fbcm = np.inner(fbsm,uBmat.T)
1633#form e.F dot product
1634            eDotFa = np.sum(eM[:,nxs,:]*facm,axis=-1)    #Nref,Ntau       
1635            eDotFb = np.sum(eM[:,nxs,:]*fbcm,axis=-1)
1636#intensity Halpern & Johnson
1637#            fass = np.sum((facm-eM[:,nxs,:]*eDotFa[:,:,nxs])**2,axis=-1)
1638#            fbss = np.sum((fbcm-eM[:,nxs,:]*eDotFb[:,:,nxs])**2,axis=-1)
1639            fass = np.sum(facm**2,axis=-1)-eDotFa**2
1640            fbss = np.sum(fbcm**2,axis=-1)-eDotFb**2
1641## #do integration           
1642            fas = np.sum(fass*glWt[nxs,:],axis=1)
1643            fbs = np.sum(fbss*glWt[nxs,:],axis=1)
1644           
1645            refl.T[10] = fas+fbs   #Sum(fams**2,Mxyz) Re + Im
1646            refl.T[11] = atan2d(fbs,fas)
1647
1648        else:
1649            GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1650            if 'T' in calcControls[hfx+'histType']:
1651                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1652                fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1653            else:
1654                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1655                fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1656            fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1657            fbg = fb*GfpuA[0]+fa*GfpuA[1]
1658            fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1659            fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1660           
1661            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1662            refl.T[11] = atan2d(fbs[0],fas[0])  #use only tau=0; ignore f' & f"
1663        if 'P' not in calcControls[hfx+'histType']:
1664            refl.T[8] = np.copy(refl.T[10])               
1665        iBeg += blkSize
1666#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
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)] = 2.*dFdua.T[3][i]
1993        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1994        dFdvDict[pfx+'AU23:'+str(i)] = 2.*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)] = 2.*dFdGu.T[3][j][i]
2013            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2014            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*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)] = 2.*dFdGu.T[9][j][i]
2019            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2020            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*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)] = 2.*dFdua.T[3][i]
2227        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2228        dFdvDict[pfx+'AU23:'+str(i)] = 2.*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)] = 2.*dFdGu.T[3][j][i]
2252            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2253            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*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)] = 2.*dFdGu.T[9][j][i]
2258            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2259            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*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 = parmDict[hfx+'alpha-0']+parmDict[hfx+'alpha-0']*tanPos
3197        bet = 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    fixedBkg = {i:Histogram['Background'][1][i] for i in Histogram['Background'][1] if i.startswith("_")} 
3204    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x,fixedBkg)
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 someting if Dij << 0.
3243        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3244        if np.any(np.diag(G)<0.) or np.any(np.isnan(A)):
3245            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
3246        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3247        Vst = np.sqrt(nl.det(G))    #V*
3248        if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
3249            if im:
3250                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
3251            elif parmDict[pfx+'isMag'] and 'N' in calcControls[hfx+'histType']:
3252                MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3253            else:
3254                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3255        badPeak = False
3256        # test to see if we are using multiprocessing here
3257        useMP,ncores = G2mp.InitMP()
3258        if len(refDict['RefList']) < 100: useMP = False       
3259        if useMP: # multiprocessing: create a set of initialized Python processes
3260            MPpool = mp.Pool(ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
3261            profArgs = [[] for i in range(ncores)]
3262        if 'C' in calcControls[hfx+'histType']:
3263            for iref,refl in enumerate(refDict['RefList']):
3264                if im:
3265                    h,k,l,m = refl[:4]
3266                else:
3267                    h,k,l = refl[:3]
3268                Uniq = np.inner(refl[:3],SGMT)
3269                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,phfx,calcControls,parmDict)         #corrected reflection position
3270                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
3271                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3272                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3273                refl[11+im] *= Vst*Lorenz
3274                 
3275                if Phase['General'].get('doPawley'):
3276                    try:
3277                        if im:
3278                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3279                        else:
3280                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3281                        refl[9+im] = parmDict[pInd]
3282                    except KeyError:
3283#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3284                        continue
3285                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
3286                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3287                iFin = np.searchsorted(x,refl[5+im]+fmax)
3288                if not iBeg+iFin:       #peak below low limit - skip peak
3289                    continue
3290                elif not iBeg-iFin:     #peak above high limit - done
3291                    break
3292                elif iBeg > iFin:   #bad peak coeff - skip
3293                    badPeak = True
3294                    continue
3295                if useMP:
3296                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin,1.))
3297                else:
3298                    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
3299                if Ka2:
3300                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3301                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
3302                    iBeg = np.searchsorted(x,pos2-fmin)
3303                    iFin = np.searchsorted(x,pos2+fmax)
3304                    if not iBeg+iFin:       #peak below low limit - skip peak
3305                        continue
3306                    elif not iBeg-iFin:     #peak above high limit - done
3307                        return yc,yb
3308                    elif iBeg > iFin:   #bad peak coeff - skip
3309                        continue
3310                    if useMP:
3311                        profArgs[iref%ncores].append((pos2,refl,iBeg,iFin,kRatio))
3312                    else:
3313                        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
3314        elif 'B' in calcControls[hfx+'histType']:
3315            for iref,refl in enumerate(refDict['RefList']):
3316                if im:
3317                    h,k,l,m = refl[:4]
3318                else:
3319                    h,k,l = refl[:3]
3320                Uniq = np.inner(refl[:3],SGMT)
3321                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,phfx,calcControls,parmDict)         #corrected reflection position
3322                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
3323                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3324                refl[12+im:14+im] = GetPinkReflAlpBet(refl,im,hfx,parmDict)             #TODO - skip if alp, bet tabulated?
3325                refl[11+im],refl[14+im],refl[15+im],refl[16+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3326                refl[11+im] *= Vst*Lorenz
3327                 
3328                if Phase['General'].get('doPawley'):
3329                    try:
3330                        if im:
3331                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3332                        else:
3333                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3334                        refl[9+im] = parmDict[pInd]
3335                    except KeyError:
3336#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3337