source: trunk/GSASIIstrMath.py @ 3372

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

disable printing of debug info for binary location & config controls
fix bug in image plotting wrt linescan
better (but not best) mag moment derivatives

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