source: trunk/GSASIIstrMath.py @ 4169

Last change on this file since 4169 was 4169, checked in by vondreele, 2 years ago

mag st fr

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