source: trunk/GSASIIstrMath.py @ 3414

Last change on this file since 3414 was 3414, checked in by vondreele, 5 years ago

correct Uij & sin/cosUij derivatives for i!=j; all needed *2.

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