source: trunk/GSASIIstrMath.py @ 3880

Last change on this file since 3880 was 3880, checked in by vondreele, 3 years ago

mag str fctr. stuff

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