source: trunk/GSASIIstrMath.py @ 3224

Last change on this file since 3224 was 3224, checked in by vondreele, 4 years ago

fix phase offset in modulations
mark places to fix mag structure factor calc & derivs.
fix mag phase import from GSAS EXP files

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