source: trunk/GSASIIstrMath.py @ 4136

Last change on this file since 4136 was 4136, checked in by toby, 3 years ago

fix cancel in refine dialog, tested w/wx2.8,3.0 & 4.0

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