source: trunk/GSASIIstrMath.py @ 4140

Last change on this file since 4140 was 4140, checked in by vondreele, 3 years ago

save some mSF work

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 219.7 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2019-09-09 16:50:58 +0000 (Mon, 09 Sep 2019) $
8# $Author: vondreele $
9# $Revision: 4140 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 4140 2019-09-09 16:50:58Z vondreele $
12########### SVN repository information ###################
13from __future__ import division, print_function
14import time
15import copy
16import numpy as np
17import numpy.ma as ma
18import numpy.linalg as nl
19import scipy.stats as st
20import multiprocessing as mp
21import GSASIIpath
22GSASIIpath.SetVersionNumber("$Revision: 4140 $")
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    g = nl.inv(G)
1485    ainv = np.sqrt(np.diag(g))
1486    GS = G/np.outer(ast,ast)
1487    uAmat = G2lat.Gmat2AB(GS)[0]
1488    Ginv = g/np.outer(ainv,ainv)
1489    Mast = twopisq*np.multiply.outer(ast,ast)   
1490    SGInv = SGData['SGInv']
1491    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1492#    Ncen = len(SGData['SGCen'])
1493    Nops = len(SGMT)*(1+SGData['SGInv'])
1494    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1495    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1496    SSCen = SSGData['SSGCen']
1497    FFtables = calcControls['FFtables']
1498    BLtables = calcControls['BLtables']
1499    MFtables = calcControls['MFtables']
1500    Amat,Bmat = G2lat.Gmat2AB(G)
1501    Flack = 1.0
1502    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1503        Flack = 1.-2.*parmDict[phfx+'Flack']
1504    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1505        GetAtomFXU(pfx,calcControls,parmDict)
1506    if not Xdata.size:          #no atoms in phase!
1507        return
1508    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1509    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)      #NB: Mmod is ReIm,Mxyz,Ntau,Natm
1510    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1511
1512    if parmDict[pfx+'isMag']:       #This part correct for making modulated mag moments on equiv atoms
1513       
1514        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
1515        MmodA,MmodB = G2mth.MagMod(glTau,mXYZ,modQ,MSSdata,SGData,SSGData)  #Ntau,Nops,Natm,Mxyz sum matches drawing
1516        Mmod = MmodA+MmodB
1517       
1518        if not SGData['SGGray']:    #for fixed Mx,My,Mz
1519            GSdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2))  #apply sym. ops.--> Natm,Nops,Nxyz
1520            if SGData['SGInv'] and not SGData['SGFixed']:   #inversion if any
1521                GSdata = np.hstack((GSdata,-GSdata))     
1522            GSdata = np.hstack([GSdata for cen in SSCen])        #dup over cell centering - Natm,Nops,Mxyz
1523            GSdata = SGData['MagMom'][nxs,:,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
1524            GSdata = np.swapaxes(GSdata,0,1)    #Nop,Natm,Mxyz
1525            Mmod += GSdata[nxs,:,:,:]
1526           
1527        Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Mmod])     #Ntau,Nop,Natm
1528        Kdata = np.inner(Mmod,uAmat)     #Cartesian vectors; Mxyz,Natm,Nop,Ntau
1529        Kmean = np.sqrt(np.sum(Kdata**2,axis=-1))
1530        Kdata /= Kmean[:,:,:,nxs]      #normalized cartesian mag. moment Ntau,Nop,Natm,Mxyz
1531
1532    FF = np.zeros(len(Tdata))
1533    if 'NC' in calcControls[hfx+'histType']:
1534        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1535    elif 'X' in calcControls[hfx+'histType']:
1536        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1537        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1538    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1539    bij = Mast*Uij
1540    blkSize = 48       #no. of reflections in a block
1541    nRef = refDict['RefList'].shape[0]
1542    SQ = 1./(2.*refDict['RefList'].T[5])**2
1543    if 'N' in calcControls[hfx+'histType']:
1544        dat = G2el.getBLvalues(BLtables)
1545        refDict['FF']['El'] = list(dat.keys())
1546        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1547        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1548        for iel,El in enumerate(refDict['FF']['El']):
1549            if El in MFtables:
1550                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1551    else:
1552        dat = G2el.getFFvalues(FFtables,0.)
1553        refDict['FF']['El'] = list(dat.keys())
1554        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1555        for iel,El in enumerate(refDict['FF']['El']):
1556            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1557#    time0 = time.time()
1558#reflection processing begins here - big arrays!
1559    iBeg = 0
1560    while iBeg < nRef:
1561        iFin = min(iBeg+blkSize,nRef)
1562        mRef= iFin-iBeg
1563        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1564        H = refl.T[:4]                          #array(blkSize,4)
1565        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1566        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1567        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1568        Uniq = np.inner(H.T,SSGMT)
1569        UniqP = np.inner(HP.T,SGMT)
1570        Phi = np.inner(H.T,SSGT)
1571        if SGInv and not SGData['SGFixed']:   #if centro - expand HKL sets
1572            Uniq = np.hstack((Uniq,-Uniq))
1573            Phi = np.hstack((Phi,-Phi))
1574            UniqP = np.hstack((UniqP,-UniqP))
1575        if 'T' in calcControls[hfx+'histType']:
1576            if 'P' in calcControls[hfx+'histType']:
1577                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1578            else:
1579                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1580            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
1581            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
1582        Bab = 0.
1583        if phfx+'BabA' in parmDict:
1584            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
1585        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1586        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
1587        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1588        phase = np.hstack([phase for cen in SSCen])
1589        sinp = np.sin(phase)
1590        cosp = np.cos(phase)
1591        biso = -SQfactor*Uisodata[:,nxs]
1592        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T
1593        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1594        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1595        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1596
1597        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:       #TODO: mag math here??           
1598           
1599            phasem = twopi*np.inner(H.T[:,:3],mXYZ)
1600            phasem = np.swapaxes(phasem,1,2)            #Nref,Nops,Natm
1601            cosm = np.cos(phasem)
1602            sinm = np.sin(phasem)
1603            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1604            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Mdata*Fdata*MF/(2.*Nops)     #Nref,Natm
1605#                     
1606            HM = np.inner(Bmat,HP.T)                            #put into cartesian space
1607            eM = (HM/np.sqrt(np.sum(HM**2,axis=0))).T               #& normalize; Nref,Mhkl
1608            eDotK = np.sum(eM[:,nxs,nxs,nxs,:]*Kdata[nxs,:,:,:,:],axis=-1)  #Nref,Ntau,Nops,Natm
1609            Q = eM[:,nxs,nxs,nxs,:]*eDotK[:,:,:,:,nxs]-Kdata[nxs,:,:,:,:]   #Nref,Ntau,Nops,Natm,Mxyz
1610            Q = 1.0
1611           
1612            fam0 = 0.
1613            fbm0 = 0.
1614            if not SGData['SGGray']:     #correct -fixed Mx,My,Mz contribution             
1615                fam0 = Q*TMcorr[:,nxs,nxs,:,nxs]*GSdata[nxs,nxs,:,:,:]*cosm[:,nxs,:,:,nxs]    #Nref,Ntau,Nops,Natm,Mxyz
1616                fbm0 = Q*TMcorr[:,nxs,nxs,:,nxs]*GSdata[nxs,nxs,:,:,:]*sinm[:,nxs,:,:,nxs]   
1617#for modulated moments --> m != 0 reflections
1618       
1619            fams = Q*TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*cosm[i,nxs,:,:,nxs]-    \
1620                np.sign(H[3,i])*MmodB*sinm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
1621                       
1622            fbms = Q*TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*sinm[i,nxs,:,:,nxs]+    \
1623                np.sign(H[3,i])*MmodB*cosm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
1624           
1625            if not SGData['SGGray']:
1626                fams *= .5
1627                fbms *= .5
1628                fams += fam0
1629                fbms += fbm0
1630               
1631#   do GL integral before making mag F's - less correct; gives F's where there shouldn't be any                     
1632#                       
1633#            fass = np.sum(glWt[nxs,:,nxs,nxs,nxs]*fams,axis=1)         #GL integration
1634#            fbss = np.sum(glWt[nxs,:,nxs,nxs,nxs]*fbms,axis=1)
1635#           
1636#            fas = np.sum(np.sum(fass,axis=-1),axis=-1)      #Nref,Ntau,Mxyz; sum ops & atoms
1637#            fbs = np.sum(np.sum(fbss,axis=-1),axis=-1)
1638#            refl.T[10] = np.sum(fas**2,axis=1)+np.sum(fbs**2,axis=1)
1639#            refl.T[11] = atan2d(fbs[:,0],fas[:,0])
1640##
1641# do GL integral after making mag. F's - better but not right               
1642               
1643            famqs = np.sum(np.sum(fams,axis=-2),axis=-2)      #Nref,Ntau,Mxyz; sum ops & atoms
1644            fbmqs = np.sum(np.sum(fbms,axis=-2),axis=-2)
1645           
1646            famcs = np.swapaxes(np.inner(Amat,famqs).T,0,1)     #convert to cartesian
1647            fbmcs = np.swapaxes(np.inner(Amat,fbmqs).T,0,1)     # as Nref,Ntau,Mxyz
1648           
1649            famcs /= np.sqrt(np.sum(famcs**2,axis=-1))[:,:,nxs]     #normalize
1650            fbmcs /= np.sqrt(np.sum(fbmcs**2,axis=-1))[:,:,nxs]
1651           
1652            famcs = np.nan_to_num(famcs)    #nan --> 0.0
1653            fbmcs = np.nan_to_num(fbmcs)
1654           
1655            fass = np.sum(famqs**2,axis=-1)*(1.-np.sum(eM[:,nxs,:]*famcs,axis=-1)**2)      #mag intensity calc F^2-(e.F)^2
1656            fbss = np.sum(fbmqs**2,axis=-1)*(1.-np.sum(eM[:,nxs,:]*fbmcs,axis=-1)**2)       #WRONG - needs mag. moment not F!
1657           
1658            fas = np.sum(glWt*fass,axis=1)          #GL integration
1659            fbs = np.sum(glWt*fbss,axis=1)
1660           
1661            refl.T[10] = fas+fbs
1662            refl.T[11] = atan2d(fbs,fas)
1663#
1664        else:
1665            GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1666            if 'T' in calcControls[hfx+'histType']:
1667                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1668                fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1669            else:
1670                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1671                fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1672            fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1673            fbg = fb*GfpuA[0]+fa*GfpuA[1]
1674            fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1675            fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1676           
1677            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1678            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1679        if 'P' not in calcControls[hfx+'histType']:
1680            refl.T[8] = np.copy(refl.T[10])               
1681        iBeg += blkSize
1682#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1683    return copy.deepcopy(refDict['RefList'])
1684
1685def SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1686    '''
1687    Compute super structure factors for all h,k,l,m for phase - twins only
1688    puts the result, F^2, in each ref[8+im] in refList
1689    works on blocks of 32 reflections for speed
1690    input:
1691   
1692    :param dict refDict: where
1693        'RefList' list where each ref = h,k,l,m,it,d,...
1694        'FF' dict of form factors - filed in below
1695    :param np.array G:      reciprocal metric tensor
1696    :param str pfx:    phase id string
1697    :param dict SGData: space group info. dictionary output from SpcGroup
1698    :param dict calcControls:
1699    :param dict ParmDict:
1700
1701    '''
1702    phfx = pfx.split(':')[0]+hfx
1703    ast = np.sqrt(np.diag(G))
1704    Mast = twopisq*np.multiply.outer(ast,ast)   
1705    SGInv = SGData['SGInv']
1706    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1707    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1708    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1709    FFtables = calcControls['FFtables']
1710    BLtables = calcControls['BLtables']
1711    MFtables = calcControls['MFtables']
1712    Flack = 1.0
1713    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1714        Flack = 1.-2.*parmDict[phfx+'Flack']
1715    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])    #4D?
1716    TwDict = refDict.get('TwDict',{})           
1717    if 'S' in calcControls[hfx+'histType']:
1718        NTL = calcControls[phfx+'NTL']
1719        NM = calcControls[phfx+'TwinNMN']+1
1720        TwinLaw = calcControls[phfx+'TwinLaw']  #this'll have to be 4D also...
1721        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1722        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1723    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1724        GetAtomFXU(pfx,calcControls,parmDict)
1725    if not Xdata.size:          #no atoms in phase!
1726        return
1727    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1728    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1729    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1730    FF = np.zeros(len(Tdata))
1731    if 'NC' in calcControls[hfx+'histType']:
1732        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1733    elif 'X' in calcControls[hfx+'histType']:
1734        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1735        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1736    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1737    bij = Mast*Uij
1738    blkSize = 32       #no. of reflections in a block
1739    nRef = refDict['RefList'].shape[0]
1740    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
1741        SQ = 1./(2.*refDict['RefList'].T[5])**2
1742        if 'N' in calcControls[hfx+'histType']:
1743            dat = G2el.getBLvalues(BLtables)
1744            refDict['FF']['El'] = list(dat.keys())
1745            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1746            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1747            for iel,El in enumerate(refDict['FF']['El']):
1748                if El in MFtables:
1749                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1750        else:
1751            dat = G2el.getFFvalues(FFtables,0.)
1752            refDict['FF']['El'] = list(dat.keys())
1753            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1754            for iel,El in enumerate(refDict['FF']['El']):
1755                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1756#    time0 = time.time()
1757#reflection processing begins here - big arrays!
1758    iBeg = 0
1759    while iBeg < nRef:
1760        iFin = min(iBeg+blkSize,nRef)
1761        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1762        H = refl[:,:4]                          #array(blkSize,4)
1763        H3 = refl[:,:3]
1764        HP = H[:,:3]+modQ[nxs,:]*H[:,3:]        #projected hklm to hkl
1765        HP = np.inner(HP,TwinLaw)             #array(blkSize,nTwins,4)
1766        H3 = np.inner(H3,TwinLaw)       
1767        TwMask = np.any(HP,axis=-1)
1768        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1769            for ir in range(blkSize):
1770                iref = ir+iBeg
1771                if iref in TwDict:
1772                    for i in TwDict[iref]:
1773                        for n in range(NTL):
1774                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1775                            H3[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1776            TwMask = np.any(HP,axis=-1)
1777        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1778        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1779        Uniq = np.inner(H,SSGMT)
1780        Uniq3 = np.inner(H3,SGMT)
1781        UniqP = np.inner(HP,SGMT)
1782        Phi = np.inner(H,SSGT)
1783        if SGInv:   #if centro - expand HKL sets
1784            Uniq = np.hstack((Uniq,-Uniq))
1785            Uniq3 = np.hstack((Uniq3,-Uniq3))
1786            Phi = np.hstack((Phi,-Phi))
1787            UniqP = np.hstack((UniqP,-UniqP))
1788        if 'T' in calcControls[hfx+'histType']:
1789            if 'P' in calcControls[hfx+'histType']:
1790                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1791            else:
1792                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1793            FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1794            FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1795        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw))
1796        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1797        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1798        phase = twopi*(np.inner(Uniq3,(dXdata.T+Xdata.T))-Phi[:,nxs,:,nxs])
1799        sinp = np.sin(phase)
1800        cosp = np.cos(phase)
1801        biso = -SQfactor*Uisodata[:,nxs]
1802        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1803        HbH = -np.sum(UniqP[:,:,:,nxs]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1804        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1805        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1806        if 'T' in calcControls[hfx+'histType']:
1807            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1808            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1809        else:
1810            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1811            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1812        GfpuA = G2mth.ModulationTw(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1813        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1814        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1815        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1816        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1817        refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
1818        refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
1819            np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
1820        refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
1821        iBeg += blkSize
1822#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1823
1824def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1825    '''
1826    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
1827    Only Fourier component are done analytically here
1828    input:
1829   
1830    :param dict refDict: where
1831        'RefList' list where each ref = h,k,l,m,it,d,...
1832        'FF' dict of form factors - filled in below
1833    :param int im: = 1 (could be eliminated)
1834    :param np.array G:      reciprocal metric tensor
1835    :param str hfx:    histogram id string
1836    :param str pfx:    phase id string
1837    :param dict SGData: space group info. dictionary output from SpcGroup
1838    :param dict SSGData: super space group info.
1839    :param dict calcControls:
1840    :param dict ParmDict:
1841   
1842    :returns: dict dFdvDict: dictionary of derivatives
1843    '''
1844    phfx = pfx.split(':')[0]+hfx
1845    ast = np.sqrt(np.diag(G))
1846    Mast = twopisq*np.multiply.outer(ast,ast)
1847    SGInv = SGData['SGInv']
1848    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1849    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1850    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1851    FFtables = calcControls['FFtables']
1852    BLtables = calcControls['BLtables']
1853    nRef = len(refDict['RefList'])
1854    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1855        GetAtomFXU(pfx,calcControls,parmDict)
1856    if not Xdata.size:          #no atoms in phase!
1857        return {}
1858    mSize = len(Mdata)  #no. atoms
1859    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1860    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
1861    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
1862    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1863    FF = np.zeros(len(Tdata))
1864    if 'NC' in calcControls[hfx+'histType']:
1865        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1866    elif 'X' in calcControls[hfx+'histType']:
1867        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1868        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1869    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1870    bij = Mast*Uij
1871    if not len(refDict['FF']):
1872        if 'N' in calcControls[hfx+'histType']:
1873            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1874        else:
1875            dat = G2el.getFFvalues(FFtables,0.)       
1876        refDict['FF']['El'] = list(dat.keys())
1877        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1878    dFdvDict = {}
1879    dFdfr = np.zeros((nRef,mSize))
1880    dFdx = np.zeros((nRef,mSize,3))
1881    dFdui = np.zeros((nRef,mSize))
1882    dFdua = np.zeros((nRef,mSize,6))
1883    dFdbab = np.zeros((nRef,2))
1884    dFdfl = np.zeros((nRef))
1885    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1886    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1887    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1888    Flack = 1.0
1889    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1890        Flack = 1.-2.*parmDict[phfx+'Flack']
1891    time0 = time.time()
1892    nRef = len(refDict['RefList'])/100
1893    for iref,refl in enumerate(refDict['RefList']):
1894        if 'T' in calcControls[hfx+'histType']:
1895            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1896        H = np.array(refl[:4])
1897        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
1898        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1899        SQfactor = 8.0*SQ*np.pi**2
1900        Bab = 0.0
1901        if phfx+'BabA' in parmDict:
1902            dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1903            Bab = parmDict[phfx+'BabA']*dBabdA
1904        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1905        FF = refDict['FF']['FF'][iref].T[Tindx]
1906        Uniq = np.inner(H,SSGMT)
1907        Phi = np.inner(H,SSGT)
1908        UniqP = np.inner(HP,SGMT)
1909        if SGInv:   #if centro - expand HKL sets
1910            Uniq = np.vstack((Uniq,-Uniq))
1911            Phi = np.hstack((Phi,-Phi))
1912            UniqP = np.vstack((UniqP,-UniqP))
1913        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1914        sinp = np.sin(phase)
1915        cosp = np.cos(phase)
1916        occ = Mdata*Fdata/Uniq.shape[0]
1917        biso = -SQfactor*Uisodata[:,nxs]
1918        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0],axis=1).T    #ops x atoms
1919        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
1920        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1921        Hij = np.array([G2lat.UijtoU6(uij) for uij in Hij])                     #atoms x 6
1922        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1923        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1924        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
1925        fotp = FPP*Tcorr            #ops x atoms
1926        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1927        dGdf,dGdx,dGdu = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1928        # GfpuA is 2 x ops x atoms
1929        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1930        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
1931        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1932        fag = fa*GfpuA[0]-fb*GfpuA[1]
1933        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1934       
1935        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
1936        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1937        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
1938        fbx = np.array([fot*cosp,-fotp*sinp])
1939        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1940        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1941        #sum below is over Uniq
1942        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
1943        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
1944        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1945        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1946        dfadui = np.sum(-SQfactor*fag,axis=1)
1947        dfbdui = np.sum(-SQfactor*fbg,axis=1)
1948        dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
1949        dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
1950        dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
1951        dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
1952        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1953        dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1954        dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1955        dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1956        dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1957        dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
1958        dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
1959        if not SGData['SGInv']:   #Flack derivative
1960            dfadfl = np.sum(-FPP*Tcorr*sinp)
1961            dfbdfl = np.sum(FPP*Tcorr*cosp)
1962        else:
1963            dfadfl = 1.0
1964            dfbdfl = 1.0
1965        SA = fas[0]+fas[1]      #float = A+A'
1966        SB = fbs[0]+fbs[1]      #float = B+B'
1967        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
1968            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1969            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
1970                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1971            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
1972                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1973            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
1974                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1975            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
1976                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1977            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
1978                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
1979            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
1980                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
1981            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
1982                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
1983        else:                       #OK, I think
1984            dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
1985            dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
1986            dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
1987            dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
1988            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1989                           
1990            dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
1991            dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
1992            dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
1993        if phfx+'BabA' in parmDict:
1994            dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1995                2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1996        #loop over atoms - each dict entry is list of derivatives for all the reflections
1997        if not iref%100 :
1998            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
1999    for i in range(len(Mdata)):     #loop over atoms
2000        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2001        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2002        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2003        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2004        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2005        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2006        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2007        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2008        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2009        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2010        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2011        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2012            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2013            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2014        nx = 0
2015        if waveTypes[i] in ['Block','ZigZag']:
2016            nx = 1 
2017        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2018            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2019            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2020            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2021            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2022            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2023            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2024        for j in range(USSdata.shape[1]):       #loop over waves
2025            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2026            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2027            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2028            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
2029            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2030            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
2031            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2032            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2033            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2034            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
2035            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2036            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2037           
2038    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
2039    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2040    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2041    return dFdvDict
2042
2043def SStructureFactorDerv2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2044    '''
2045    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
2046    input:
2047   
2048    :param dict refDict: where
2049        'RefList' list where each ref = h,k,l,m,it,d,...
2050        'FF' dict of form factors - filled in below
2051    :param int im: = 1 (could be eliminated)
2052    :param np.array G:      reciprocal metric tensor
2053    :param str hfx:    histogram id string
2054    :param str pfx:    phase id string
2055    :param dict SGData: space group info. dictionary output from SpcGroup
2056    :param dict SSGData: super space group info.
2057    :param dict calcControls:
2058    :param dict ParmDict:
2059   
2060    :returns: dict dFdvDict: dictionary of derivatives
2061    '''
2062
2063    trefDict = copy.deepcopy(refDict)
2064    dM = 1.e-4
2065    dFdvDict = {}
2066    for parm in parmDict:
2067        if ':' not in parm:
2068            continue
2069        if parm.split(':')[2] in ['Tmin','Tmax','Xmax','Ymax','Zmax','Fzero','Fwid',
2070            'MXsin','MXcos','MYsin','MYcos','MZsin','MZcos','AMx','AMy','AMz',]:
2071            parmDict[parm] += dM
2072            prefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2073            parmDict[parm] -= 2*dM
2074            mrefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2075            parmDict[parm] += dM
2076            dFdvDict[parm] = (prefList[:,9+im]-mrefList[:,9+im])/(2.*dM)
2077    return dFdvDict
2078   
2079def SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2080    'Needs a doc string'
2081    phfx = pfx.split(':')[0]+hfx
2082    ast = np.sqrt(np.diag(G))
2083    Mast = twopisq*np.multiply.outer(ast,ast)
2084    SGInv = SGData['SGInv']
2085    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2086    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2087    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2088    FFtables = calcControls['FFtables']
2089    BLtables = calcControls['BLtables']
2090    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
2091    TwDict = refDict.get('TwDict',{})           
2092    if 'S' in calcControls[hfx+'histType']:
2093        NTL = calcControls[phfx+'NTL']
2094        NM = calcControls[phfx+'TwinNMN']+1
2095        TwinLaw = calcControls[phfx+'TwinLaw']
2096        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
2097    nTwin = len(TwinLaw)       
2098    nRef = len(refDict['RefList'])
2099    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
2100        GetAtomFXU(pfx,calcControls,parmDict)
2101    if not Xdata.size:          #no atoms in phase!
2102        return {}
2103    mSize = len(Mdata)  #no. atoms
2104    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
2105    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)     #NB: Mmod is ReIm,Mxyz,Ntau,Natm
2106    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
2107    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2108    FF = np.zeros(len(Tdata))
2109    if 'NC' in calcControls[hfx+'histType']:
2110        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
2111    elif 'X' in calcControls[hfx+'histType']:
2112        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2113        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2114    Uij = np.array(G2lat.U6toUij(Uijdata)).T
2115    bij = Mast*Uij
2116    if not len(refDict['FF']):
2117        if 'N' in calcControls[hfx+'histType']:
2118            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
2119        else:
2120            dat = G2el.getFFvalues(FFtables,0.)       
2121        refDict['FF']['El'] = list(dat.keys())
2122        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
2123    dFdvDict = {}
2124    dFdfr = np.zeros((nRef,nTwin,mSize))
2125    dFdx = np.zeros((nRef,nTwin,mSize,3))
2126    dFdui = np.zeros((nRef,nTwin,mSize))
2127    dFdua = np.zeros((nRef,nTwin,mSize,6))
2128    dFdbab = np.zeros((nRef,nTwin,2))
2129    dFdtw = np.zeros((nRef,nTwin))
2130    dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
2131    dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
2132    dFdGz = np.zeros((nRef,nTwin,mSize,5))
2133    dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
2134    Flack = 1.0
2135    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
2136        Flack = 1.-2.*parmDict[phfx+'Flack']
2137    time0 = time.time()
2138    nRef = len(refDict['RefList'])/100
2139    for iref,refl in enumerate(refDict['RefList']):
2140        if 'T' in calcControls[hfx+'histType']:
2141            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
2142        H = np.array(refl[:4])
2143        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
2144        H = np.inner(H.T,TwinLaw)   #maybe array(4,nTwins) or (4)
2145        TwMask = np.any(H,axis=-1)
2146        if TwinLaw.shape[0] > 1 and TwDict:
2147            if iref in TwDict:
2148                for i in TwDict[iref]:
2149                    for n in range(NTL):
2150                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
2151            TwMask = np.any(H,axis=-1)
2152        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
2153        SQfactor = 8.0*SQ*np.pi**2
2154        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2155        Bab = parmDict[phfx+'BabA']*dBabdA
2156        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
2157        FF = refDict['FF']['FF'][iref].T[Tindx]
2158        Uniq = np.inner(H,SSGMT)
2159        Phi = np.inner(H,SSGT)
2160        UniqP = np.inner(HP,SGMT)
2161        if SGInv:   #if centro - expand HKL sets
2162            Uniq = np.vstack((Uniq,-Uniq))
2163            Phi = np.hstack((Phi,-Phi))
2164            UniqP = np.vstack((UniqP,-UniqP))
2165        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
2166        sinp = np.sin(phase)
2167        cosp = np.cos(phase)
2168        occ = Mdata*Fdata/Uniq.shape[0]
2169        biso = -SQfactor*Uisodata[:,nxs]
2170        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
2171        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
2172        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
2173        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(nTwin,-1,6)))
2174        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
2175        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
2176        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
2177        fotp = FPP*Tcorr            #ops x atoms
2178        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
2179        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
2180        # GfpuA is 2 x ops x atoms
2181        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
2182        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
2183        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
2184        fag = fa*GfpuA[0]-fb*GfpuA[1]
2185        fbg = fb*GfpuA[0]+fa*GfpuA[1]
2186       
2187        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
2188        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
2189        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x twin x ops x atoms
2190        fbx = np.array([fot*cosp,-fotp*sinp])
2191        fax = fax*GfpuA[0]-fbx*GfpuA[1]
2192        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
2193        #sum below is over Uniq
2194        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
2195        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
2196        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
2197        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
2198        dfadui = np.sum(-SQfactor*fag,axis=1)
2199        dfbdui = np.sum(-SQfactor*fbg,axis=1)
2200        dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2201        dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
2202        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2203        dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2204        # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
2205        dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2206        dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2207        dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2208        dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2209        dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2210        dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2211        dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2212        dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2213#        GSASIIpath.IPyBreak()
2214        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
2215        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
2216        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
2217        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)]
2218        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)]
2219        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)]
2220        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)]
2221        dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
2222
2223        dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
2224        dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
2225        dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
2226        dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
2227#            GSASIIpath.IPyBreak()
2228        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
2229            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2230        #loop over atoms - each dict entry is list of derivatives for all the reflections
2231        if not iref%100 :
2232            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
2233    for i in range(len(Mdata)):     #loop over atoms
2234        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2235        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2236        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2237        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2238        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2239        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2240        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2241        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2242        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2243        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2244        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2245        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2246            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2247            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2248        nx = 0
2249        if waveTypes[i] in ['Block','ZigZag']:
2250            nx = 1 
2251            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
2252            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
2253            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
2254            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
2255            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
2256        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2257            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2258            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2259            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2260            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2261            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2262            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2263        for j in range(USSdata.shape[1]):       #loop over waves
2264            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2265            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2266            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2267            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
2268            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2269            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
2270            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2271            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2272            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2273            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
2274            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2275            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2276           
2277#        GSASIIpath.IPyBreak()
2278    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2279    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2280    return dFdvDict
2281   
2282def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
2283    ''' Single crystal extinction function; returns extinction & derivative
2284    '''
2285    extCor = 1.0
2286    dervDict = {}
2287    dervCor = 1.0
2288    if calcControls[phfx+'EType'] != 'None':
2289        SQ = 1/(4.*ref[4+im]**2)
2290        if 'C' in parmDict[hfx+'Type']:           
2291            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
2292        else:   #'T'
2293            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
2294        if 'SXC' in parmDict[hfx+'Type']:
2295            AV = 7.9406e5/parmDict[pfx+'Vol']**2
2296            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2297            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
2298            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
2299        elif 'SNT' in parmDict[hfx+'Type']:
2300            AV = 1.e7/parmDict[pfx+'Vol']**2
2301            PL = SQ
2302            PLZ = AV*ref[9+im]*ref[12+im]**2
2303        elif 'SNC' in parmDict[hfx+'Type']:
2304            AV = 1.e7/parmDict[pfx+'Vol']**2
2305            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2306            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
2307           
2308        if 'Primary' in calcControls[phfx+'EType']:
2309            PLZ *= 1.5
2310        else:
2311            if 'C' in parmDict[hfx+'Type']:
2312                PLZ *= calcControls[phfx+'Tbar']
2313            else: #'T'
2314                PLZ *= ref[13+im]      #t-bar
2315        if 'Primary' in calcControls[phfx+'EType']:
2316            PLZ *= 1.5
2317            PSIG = parmDict[phfx+'Ep']
2318        elif 'I & II' in calcControls[phfx+'EType']:
2319            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
2320        elif 'Type II' in calcControls[phfx+'EType']:
2321            PSIG = parmDict[phfx+'Es']
2322        else:       # 'Secondary Type I'
2323            PSIG = parmDict[phfx+'Eg']/PL
2324           
2325        AG = 0.58+0.48*cos2T+0.24*cos2T**2
2326        AL = 0.025+0.285*cos2T
2327        BG = 0.02-0.025*cos2T
2328        BL = 0.15-0.2*(0.75-cos2T)**2
2329        if cos2T < 0.:
2330            BL = -0.45*cos2T
2331        CG = 2.
2332        CL = 2.
2333        PF = PLZ*PSIG
2334       
2335        if 'Gaussian' in calcControls[phfx+'EApprox']:
2336            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
2337            extCor = np.sqrt(PF4)
2338            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
2339        else:
2340            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
2341            extCor = np.sqrt(PF4)
2342            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
2343
2344        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
2345        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
2346            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
2347        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
2348            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
2349        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
2350            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
2351               
2352    return 1./extCor,dervDict,dervCor
2353   
2354def Dict2Values(parmdict, varylist):
2355    '''Use before call to leastsq to setup list of values for the parameters
2356    in parmdict, as selected by key in varylist'''
2357    return [parmdict[key] for key in varylist] 
2358   
2359def Values2Dict(parmdict, varylist, values):
2360    ''' Use after call to leastsq to update the parameter dictionary with
2361    values corresponding to keys in varylist'''
2362    parmdict.update(zip(varylist,values))
2363   
2364def GetNewCellParms(parmDict,varyList):
2365    '''Compute unit cell tensor terms from varied Aij and Dij values.
2366    Terms are included in the dict only if Aij or Dij is varied.
2367    '''
2368    newCellDict = {}
2369    Anames = ['A'+str(i) for i in range(6)]
2370    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
2371    for item in varyList:
2372        keys = item.split(':')
2373        if keys[2] in Ddict:
2374            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
2375            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
2376            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2377    return newCellDict          # is e.g. {'0::D11':A0-D11}
2378   
2379def ApplyXYZshifts(parmDict,varyList):
2380    '''
2381    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2382   
2383    :param dict parmDict: parameter dictionary
2384    :param list varyList: list of variables (not used!)
2385    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
2386
2387    '''
2388    newAtomDict = {}
2389    for item in parmDict:
2390        if 'dA' in item:
2391            parm = ''.join(item.split('d'))
2392            parmDict[parm] += parmDict[item]
2393            newAtomDict[item] = [parm,parmDict[parm]]
2394    return newAtomDict
2395   
2396def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2397    'Spherical harmonics texture'
2398    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2399    if 'T' in calcControls[hfx+'histType']:
2400        tth = parmDict[hfx+'2-theta']
2401    else:
2402        tth = refl[5+im]
2403    odfCor = 1.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,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
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,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2415        Lnorm = G2lat.Lnorm(L)
2416        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2417    return odfCor
2418   
2419def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2420    'Spherical harmonics texture derivatives'
2421    if 'T' in calcControls[hfx+'histType']:
2422        tth = parmDict[hfx+'2-theta']
2423    else:
2424        tth = refl[5+im]
2425    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2426    odfCor = 1.0
2427    dFdODF = {}
2428    dFdSA = [0,0,0]
2429    H = refl[:3]
2430    cell = G2lat.Gmat2cell(g)
2431    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2432    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2433    phi,beta = G2lat.CrsAng(H,cell,SGData)
2434    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
2435    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2436    for item in SHnames:
2437        L,M,N = eval(item.strip('C'))
2438        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2439        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2440        Lnorm = G2lat.Lnorm(L)
2441        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2442        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2443        for i in range(3):
2444            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2445    return odfCor,dFdODF,dFdSA
2446   
2447def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2448    'spherical harmonics preferred orientation (cylindrical symmetry only)'
2449    if 'T' in calcControls[hfx+'histType']:
2450        tth = parmDict[hfx+'2-theta']
2451    else:
2452        tth = refl[5+im]
2453    odfCor = 1.0
2454    H = refl[:3]
2455    cell = G2lat.Gmat2cell(g)
2456    Sangls = [0.,0.,0.]
2457    if 'Bragg' in calcControls[hfx+'instType']:
2458        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2459        IFCoup = True
2460    else:
2461        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2462        IFCoup = False
2463    phi,beta = G2lat.CrsAng(H,cell,SGData)
2464    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2465    SHnames = calcControls[phfx+'SHnames']
2466    for item in SHnames:
2467        L,N = eval(item.strip('C'))
2468        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2469        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2470        Lnorm = G2lat.Lnorm(L)
2471        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2472    return np.squeeze(odfCor)
2473   
2474def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2475    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
2476    if 'T' in calcControls[hfx+'histType']:
2477        tth = parmDict[hfx+'2-theta']
2478    else:
2479        tth = refl[5+im]
2480    odfCor = 1.0
2481    dFdODF = {}
2482    H = refl[:3]
2483    cell = G2lat.Gmat2cell(g)
2484    Sangls = [0.,0.,0.]
2485    if 'Bragg' in calcControls[hfx+'instType']:
2486        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2487        IFCoup = True
2488    else:
2489        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2490        IFCoup = False
2491    phi,beta = G2lat.CrsAng(H,cell,SGData)
2492    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2493    SHnames = calcControls[phfx+'SHnames']
2494    for item in SHnames:
2495        L,N = eval(item.strip('C'))
2496        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2497        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2498        Lnorm = G2lat.Lnorm(L)
2499        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2500        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
2501    return odfCor,dFdODF
2502   
2503def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2504    'March-Dollase preferred orientation correction'
2505    POcorr = 1.0
2506    MD = parmDict[phfx+'MD']
2507    if MD != 1.0:
2508        MDAxis = calcControls[phfx+'MDAxis']
2509        sumMD = 0
2510        for H in uniq:           
2511            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2512            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2513            sumMD += A**3
2514        POcorr = sumMD/len(uniq)
2515    return POcorr
2516   
2517def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2518    'Needs a doc string'
2519    POcorr = 1.0
2520    POderv = {}
2521    if calcControls[phfx+'poType'] == 'MD':
2522        MD = parmDict[phfx+'MD']
2523        MDAxis = calcControls[phfx+'MDAxis']
2524        sumMD = 0
2525        sumdMD = 0
2526        for H in uniq:           
2527            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2528            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2529            sumMD += A**3
2530            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2531        POcorr = sumMD/len(uniq)
2532        POderv[phfx+'MD'] = sumdMD/len(uniq)
2533    else:   #spherical harmonics
2534        if calcControls[phfx+'SHord']:
2535            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2536    return POcorr,POderv
2537   
2538def GetAbsorb(refl,im,hfx,calcControls,parmDict):
2539    'Needs a doc string'
2540    if 'Debye' in calcControls[hfx+'instType']:
2541        if 'T' in calcControls[hfx+'histType']:
2542            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2543        else:
2544            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2545    else:
2546        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
2547   
2548def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
2549    'Needs a doc string'
2550    if 'Debye' in calcControls[hfx+'instType']:
2551        if 'T' in calcControls[hfx+'histType']:
2552            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2553        else:
2554            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2555    else:
2556        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
2557       
2558def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2559    'Needs a doc string'
2560    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2561    pi2 = np.sqrt(2./np.pi)
2562    if 'T' in calcControls[hfx+'histType']:
2563        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2564        wave = refl[14+im]
2565    else:   #'C'W
2566        sth2 = sind(refl[5+im]/2.)**2
2567        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2568    c2th = 1.-2.0*sth2
2569    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2570    if 'X' in calcControls[hfx+'histType']:
2571        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2572    xfac = flv2*parmDict[phfx+'Extinction']
2573    exb = 1.0
2574    if xfac > -1.:
2575        exb = 1./np.sqrt(1.+xfac)
2576    exl = 1.0
2577    if 0 < xfac <= 1.:
2578        xn = np.array([xfac**(i+1) for i in range(6)])
2579        exl += np.sum(xn*coef)
2580    elif xfac > 1.:
2581        xfac2 = 1./np.sqrt(xfac)
2582        exl = pi2*(1.-0.125/xfac)*xfac2
2583    return exb*sth2+exl*(1.-sth2)
2584   
2585def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2586    'Needs a doc string'
2587    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2588    pi2 = np.sqrt(2./np.pi)
2589    if 'T' in calcControls[hfx+'histType']:
2590        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2591        wave = refl[14+im]
2592    else:   #'C'W
2593        sth2 = sind(refl[5+im]/2.)**2
2594        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2595    c2th = 1.-2.0*sth2
2596    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2597    if 'X' in calcControls[hfx+'histType']:
2598        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2599    xfac = flv2*parmDict[phfx+'Extinction']
2600    dbde = -500.*flv2
2601    if xfac > -1.:
2602        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
2603    dlde = 0.
2604    if 0 < xfac <= 1.:
2605        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
2606        dlde = np.sum(xn*coef)/xfac
2607    elif xfac > 1.:
2608        xfac2 = 1./np.sqrt(xfac)
2609        dlde = 0.5*flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
2610       
2611    return dbde*sth2+dlde*(1.-sth2)
2612   
2613def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2614    'Needs a doc string'    #need powder extinction!
2615    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
2616    if 'X' in parmDict[hfx+'Type']:
2617        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
2618    POcorr = 1.0
2619    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
2620        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2621    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
2622        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
2623    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
2624        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2625    Icorr *= POcorr
2626    AbsCorr = 1.0
2627    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
2628    Icorr *= AbsCorr
2629    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
2630    Icorr *= ExtCorr
2631    return Icorr,POcorr,AbsCorr,ExtCorr
2632   
2633def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2634    'Needs a doc string'    #need powder extinction derivs!
2635    dIdsh = 1./parmDict[hfx+'Scale']
2636    dIdsp = 1./parmDict[phfx+'Scale']
2637    if 'X' in parmDict[hfx+'Type']:
2638        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
2639        dIdPola /= pola
2640    else:       #'N'
2641        dIdPola = 0.0
2642    dFdODF = {}
2643    dFdSA = [0,0,0]
2644    dIdPO = {}
2645    if pfx+'SHorder' in parmDict:
2646        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2647        for iSH in dFdODF:
2648            dFdODF[iSH] /= odfCor
2649        for i in range(3):
2650            dFdSA[i] /= odfCor
2651    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
2652        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
2653        for iPO in dIdPO:
2654            dIdPO[iPO] /= POcorr
2655    if 'T' in parmDict[hfx+'Type']:
2656        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
2657        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
2658    else:
2659        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
2660        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
2661    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
2662       
2663def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2664    'Needs a doc string'
2665    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
2666        costh = cosd(refl[5+im]/2.)
2667        #crystallite size
2668        if calcControls[phfx+'SizeType'] == 'isotropic':
2669            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2670        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2671            H = np.array(refl[:3])
2672            P = np.array(calcControls[phfx+'SizeAxis'])
2673            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2674            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
2675            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2676        else:           #ellipsoidal crystallites
2677            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2678            H = np.array(refl[:3])
2679            lenR = G2pwd.ellipseSize(H,Sij,GB)
2680            Sgam = 1.8*wave/(np.pi*costh*lenR)
2681        #microstrain               
2682        if calcControls[phfx+'MustrainType'] == 'isotropic':
2683            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2684        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2685            H = np.array(refl[:3])
2686            P = np.array(calcControls[phfx+'MustrainAxis'])
2687            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2688            Si = parmDict[phfx+'Mustrain;i']
2689            Sa = parmDict[phfx+'Mustrain;a']
2690            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2691        else:       #generalized - P.W. Stephens model
2692            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2693            Sum = 0
2694            for i,strm in enumerate(Strms):
2695                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2696            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
2697    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
2698        #crystallite size
2699        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2700            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2701        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2702            H = np.array(refl[:3])
2703            P = np.array(calcControls[phfx+'SizeAxis'])
2704            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2705            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
2706            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2707        else:           #ellipsoidal crystallites   #OK
2708            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2709            H = np.array(refl[:3])
2710            lenR = G2pwd.ellipseSize(H,Sij,GB)
2711            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
2712        #microstrain               
2713        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
2714            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2715        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
2716            H = np.array(refl[:3])
2717            P = np.array(calcControls[phfx+'MustrainAxis'])
2718            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2719            Si = parmDict[phfx+'Mustrain;i']
2720            Sa = parmDict[phfx+'Mustrain;a']
2721            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2722        else:       #generalized - P.W. Stephens model  OK
2723            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2724            Sum = 0
2725            for i,strm in enumerate(Strms):
2726                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2727            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
2728           
2729    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
2730    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
2731    sig /= ateln2
2732    return sig,gam
2733       
2734def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2735    'Needs a doc string'
2736    gamDict = {}
2737    sigDict = {}
2738    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
2739        costh = cosd(refl[5+im]/2.)
2740        tanth = tand(refl[5+im]/2.)
2741        #crystallite size derivatives
2742        if calcControls[phfx+'SizeType'] == 'isotropic':
2743            Sgam = 1.8*wave/(np.pi*costh*parmDict[phfx+'Size;i'])
2744            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh*parmDict[phfx+'Size;i']**2)
2745            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
2746        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2747            H = np.array(refl[:3])
2748            P = np.array(calcControls[phfx+'SizeAxis'])
2749            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2750            Si = parmDict[phfx+'Size;i']
2751            Sa = parmDict[phfx+'Size;a']
2752            gami = 1.8*wave/(costh*np.pi*Si*Sa)
2753            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2754            Sgam = gami*sqtrm
2755            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2756            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2757            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2758            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2759            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2760            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2761        else:           #ellipsoidal crystallites
2762            const = 1.8*wave/(np.pi*costh)
2763            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2764            H = np.array(refl[:3])
2765            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2766            Sgam = const/lenR
2767            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2768                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2769                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2770        gamDict[phfx+'Size;mx'] = Sgam
2771        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
2772               
2773        #microstrain derivatives               
2774        if calcControls[phfx+'MustrainType'] == 'isotropic':
2775            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2776            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
2777            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
2778        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2779            H = np.array(refl[:3])
2780            P = np.array(calcControls[phfx+'MustrainAxis'])
2781            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2782            Si = parmDict[phfx+'Mustrain;i']
2783            Sa = parmDict[phfx+'Mustrain;a']
2784            gami = 0.018*Si*Sa*tanth/np.pi
2785            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2786            Mgam = gami/sqtrm
2787            dsi = -gami*Si*cosP**2/sqtrm**3
2788            dsa = -gami*Sa*sinP**2/sqtrm**3
2789            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2790            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2791            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2792            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2793        else:       #generalized - P.W. Stephens model
2794            const = 0.018*refl[4+im]**2*tanth/np.pi
2795            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2796            Sum = 0
2797            for i,strm in enumerate(Strms):
2798                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2799                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2800                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2801            Mgam = const*np.sqrt(Sum)
2802            for i in range(len(Strms)):
2803                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2804                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2
2805        gamDict[phfx+'Mustrain;mx'] = Mgam
2806        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2807    else:   #'T'OF - All checked & OK
2808        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2809            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2810            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
2811            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
2812        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2813            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2814            H = np.array(refl[:3])
2815            P = np.array(calcControls[phfx+'SizeAxis'])
2816            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2817            Si = parmDict[phfx+'Size;i']
2818            Sa = parmDict[phfx+'Size;a']
2819            gami = const/(Si*Sa)
2820            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2821            Sgam = gami*sqtrm
2822            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2823            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2824            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2825            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2826            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2827            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2828        else:           #OK  ellipsoidal crystallites
2829            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2830            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2831            H = np.array(refl[:3])
2832            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2833            Sgam = const/lenR
2834            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2835                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2836                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2837        gamDict[phfx+'Size;mx'] = Sgam  #OK
2838        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
2839               
2840        #microstrain derivatives               
2841        if calcControls[phfx+'MustrainType'] == 'isotropic':
2842            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2843            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
2844            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
2845        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2846            H = np.array(refl[:3])
2847            P = np.array(calcControls[phfx+'MustrainAxis'])
2848            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2849            Si = parmDict[phfx+'Mustrain;i']
2850            Sa = parmDict[phfx+'Mustrain;a']
2851            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
2852            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2853            Mgam = gami/sqtrm
2854            dsi = -gami*Si*cosP**2/sqtrm**3
2855            dsa = -gami*Sa*sinP**2/sqtrm**3
2856            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2857            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2858            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2859            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2860        else:       #generalized - P.W. Stephens model OK
2861            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2862            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
2863            Sum = 0
2864            for i,strm in enumerate(Strms):
2865                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2866                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2867                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2868            Mgam = const*np.sqrt(Sum)
2869            for i in range(len(Strms)):
2870                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2871                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2       
2872        gamDict[phfx+'Mustrain;mx'] = Mgam
2873        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2874       
2875    return sigDict,gamDict
2876       
2877def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2878    'Needs a doc string'
2879    if im:
2880        h,k,l,m = refl[:4]
2881        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2882        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
2883    else:
2884        h,k,l = refl[:3]
2885        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
2886    refl[4+im] = d
2887    if 'C' in calcControls[hfx+'histType']:
2888        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2889        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2890        if 'Bragg' in calcControls[hfx+'instType']:
2891            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2892                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2893        else:               #Debye-Scherrer - simple but maybe not right
2894            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2895    elif 'T' in calcControls[hfx+'histType']:
2896        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
2897        #do I need sample position effects - maybe?
2898    return pos
2899
2900def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2901    'Needs a doc string'
2902    dpr = 180./np.pi
2903    if im:
2904        h,k,l,m = refl[:4]
2905        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2906        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
2907        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
2908    else:
2909        m = 0
2910        h,k,l = refl[:3]       
2911        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2912    dst = np.sqrt(dstsq)
2913    dsp = 1./dst
2914    if 'C' in calcControls[hfx+'histType']:
2915        pos = refl[5+im]-parmDict[hfx+'Zero']
2916        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2917        dpdw = const*dst
2918        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
2919        dpdZ = 1.0
2920        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2921            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
2922        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2923        if 'Bragg' in calcControls[hfx+'instType']:
2924            dpdSh = -4.*shft*cosd(pos/2.0)
2925            dpdTr = -shft*sind(pos)*100.0
2926            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
2927        else:               #Debye-Scherrer - simple but maybe not right
2928            dpdXd = -shft*cosd(pos)
2929            dpdYd = -shft*sind(pos)
2930            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
2931    elif 'T' in calcControls[hfx+'histType']:
2932        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
2933        dpdZ = 1.0
2934        dpdDC = dsp
2935        dpdDA = dsp**2
2936        dpdDB = 1./dsp
2937        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2938            2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2.
2939        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2940           
2941def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2942    'Needs a doc string'
2943    laue = SGData['SGLaue']
2944    uniq = SGData['SGUniq']
2945    h,k,l = refl[:3]
2946    if laue in ['m3','m3m']:
2947        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2948            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2949    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2950        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2951    elif laue in ['3R','3mR']:
2952        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2953    elif laue in ['4/m','4/mmm']:
2954        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2955    elif laue in ['mmm']:
2956        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2957    elif laue in ['2/m']:
2958        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2959        if uniq == 'a':
2960            Dij += parmDict[phfx+'D23']*k*l
2961        elif uniq == 'b':
2962            Dij += parmDict[phfx+'D13']*h*l
2963        elif uniq == 'c':
2964            Dij += parmDict[phfx+'D12']*h*k
2965    else:
2966        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2967            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2968    if 'C' in calcControls[hfx+'histType']:
2969        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2970    else:
2971        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
2972           
2973def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2974    'Needs a doc string'
2975    laue = SGData['SGLaue']
2976    uniq = SGData['SGUniq']
2977    h,k,l = refl[:3]
2978    if laue in ['m3','m3m']:
2979        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2980            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2981    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2982        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2983    elif laue in ['3R','3mR']:
2984        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2985    elif laue in ['4/m','4/mmm']:
2986        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2987    elif laue in ['mmm']:
2988        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2989    elif laue in ['2/m']:
2990        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2991        if uniq == 'a':
2992            dDijDict[phfx+'D23'] = k*l
2993        elif uniq == 'b':
2994            dDijDict[phfx+'D13'] = h*l
2995        elif uniq == 'c':
2996            dDijDict[phfx+'D12'] = h*k
2997    else:
2998        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2999            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
3000    if 'C' in calcControls[hfx+'histType']:
3001        for item in dDijDict:
3002            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
3003    else:
3004        for item in dDijDict:
3005            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
3006    return dDijDict
3007   
3008def GetDij(phfx,SGData,parmDict):
3009    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
3010    return G2spc.HStrainVals(HSvals,SGData)
3011               
3012def GetFobsSq(Histograms,Phases,parmDict,calcControls):
3013    '''Compute the observed structure factors for Powder histograms and store in reflection array
3014    Multiprocessing support added
3015    '''
3016    if GSASIIpath.GetConfigValue('Show_timing',False):
3017        starttime = time.time() #; print 'start GetFobsSq'
3018    histoList = list(Histograms.keys())
3019    histoList.sort()
3020    Ka2 = shl = lamRatio = kRatio = None
3021    for histogram in histoList:
3022        if 'PWDR' in histogram[:4]:
3023            Histogram = Histograms[histogram]
3024            hId = Histogram['hId']
3025            hfx = ':%d:'%(hId)
3026            Limits = calcControls[hfx+'Limits']
3027            if 'C' in calcControls[hfx+'histType']:
3028                shl = max(parmDict[hfx+'SH/L'],0.0005)
3029                Ka2 = False
3030                kRatio = 0.0
3031                if hfx+'Lam1' in list(parmDict.keys()):
3032                    Ka2 = True
3033                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3034                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
3035            x,y,w,yc,yb,yd = Histogram['Data']
3036            xMask = ma.getmaskarray(x)
3037            xB = np.searchsorted(x,Limits[0])
3038            xF = np.searchsorted(x,Limits[1])
3039            ymb = np.array(y-yb)
3040            ymb = np.where(ymb,ymb,1.0)
3041            ycmb = np.array(yc-yb)
3042            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
3043            refLists = Histogram['Reflection Lists']
3044            for phase in refLists:
3045                if phase not in Phases:     #skips deleted or renamed phases silently!
3046                    continue
3047                Phase = Phases[phase]
3048                if histogram not in Phase['Histograms']:
3049                    continue
3050                im = 0
3051                if Phase['General'].get('Modulated',False):
3052                    im = 1
3053                pId = Phase['pId']
3054                phfx = '%d:%d:'%(pId,hId)
3055                refDict = refLists[phase]
3056                sumFo = 0.0
3057                sumdF = 0.0
3058                sumFosq = 0.0
3059                sumdFsq = 0.0
3060                sumInt = 0.0
3061                nExcl = 0
3062                # test to see if we are using multiprocessing below
3063                useMP,ncores = G2mp.InitMP()
3064                if len(refDict['RefList']) < 100: useMP = False       
3065                if useMP: # multiprocessing: create a set of initialized Python processes
3066                    MPpool = mp.Pool(G2mp.ncores,G2mp.InitFobsSqGlobals,
3067                                    [x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2])
3068                    profArgs = [[] for i in range(G2mp.ncores)]
3069                else:
3070                    G2mp.InitFobsSqGlobals(x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2)
3071                if 'C' in calcControls[hfx+'histType']:
3072                    # are we multiprocessing?
3073                    for iref,refl in enumerate(refDict['RefList']):
3074                        if useMP: 
3075                            profArgs[iref%G2mp.ncores].append((refl,iref))
3076                        else:
3077                            icod= G2mp.ComputeFobsSqCW(refl,iref)
3078                            if type(icod) is tuple:
3079                                refl[8+im] = icod[0]
3080                                sumInt += icod[1]
3081                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3082                            elif icod == -1:
3083                                refl[3+im] *= -1
3084                                nExcl += 1
3085                            elif icod == -2:
3086                                break
3087                    if useMP:
3088                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):
3089                            sumInt += sInt
3090                            for refl8im,irefl in resList:
3091                                if refl8im is None:
3092                                    refDict['RefList'][irefl][3+im] *= -1
3093                                    nExcl += 1
3094                                else:
3095                                    refDict['RefList'][irefl][8+im] = refl8im
3096                                    if parmDict[phfx+'LeBail']:
3097                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3098                elif 'T' in calcControls[hfx+'histType']:
3099                    for iref,refl in enumerate(refDict['RefList']):
3100                        if useMP: 
3101                            profArgs[iref%G2mp.ncores].append((refl,iref))
3102                        else:
3103                            icod= G2mp.ComputeFobsSqTOF(refl,iref)
3104                            if type(icod) is tuple:
3105                                refl[8+im] = icod[0]
3106                                sumInt += icod[1]
3107                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3108                            elif icod == -1:
3109                                refl[3+im] *= -1
3110                                nExcl += 1
3111                            elif icod == -2:
3112                                break
3113                    if useMP:
3114                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):
3115                            sumInt += sInt
3116                            for refl8im,irefl in resList:
3117                                if refl8im is None:
3118                                    refDict['RefList'][irefl][3+im] *= -1
3119                                    nExcl += 1
3120                                else:
3121                                    refDict['RefList'][irefl][8+im] = refl8im
3122                                    if parmDict[phfx+'LeBail']:
3123                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3124                if useMP: MPpool.terminate()
3125                sumFo = 0.0
3126                sumdF = 0.0
3127                sumFosq = 0.0
3128                sumdFsq = 0.0
3129                for iref,refl in enumerate(refDict['RefList']):
3130                    Fo = np.sqrt(np.abs(refl[8+im]))
3131                    Fc = np.sqrt(np.abs(refl[9]+im))
3132                    sumFo += Fo
3133                    sumFosq += refl[8+im]**2
3134                    sumdF += np.abs(Fo-Fc)
3135                    sumdFsq += (refl[8+im]-refl[9+im])**2
3136                if sumFo:
3137                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3138                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3139                else:
3140                    Histogram['Residuals'][phfx+'Rf'] = 100.
3141                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
3142                Histogram['Residuals'][phfx+'sumInt'] = sumInt
3143                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])-nExcl
3144                Histogram['Residuals']['hId'] = hId
3145        elif 'HKLF' in histogram[:4]:
3146            Histogram = Histograms[histogram]
3147            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
3148    if GSASIIpath.GetConfigValue('Show_timing',False):
3149        print ('GetFobsSq t=',time.time()-starttime)
3150               
3151def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup,histogram=None):
3152    'Computes the powder pattern for a histogram based on contributions from all used phases'
3153    if GSASIIpath.GetConfigValue('Show_timing',False): starttime = time.time()
3154   
3155    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
3156        U = parmDict[hfx+'U']
3157        V = parmDict[hfx+'V']
3158        W = parmDict[hfx+'W']
3159        X = parmDict[hfx+'X']
3160        Y = parmDict[hfx+'Y']
3161        Z = parmDict[hfx+'Z']
3162        tanPos = tand(refl[5+im]/2.0)
3163        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3164        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3165        sig = max(0.001,sig)
3166        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam+Z     #save peak gamma
3167        gam = max(0.001,gam)
3168        return sig,gam
3169               
3170    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
3171        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
3172            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']*refl[4+im]
3173        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2+parmDict[hfx+'Z']
3174        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3175        sig += Ssig
3176        gam += Sgam
3177        return sig,gam
3178       
3179    def GetReflAlpBet(refl,im,hfx,parmDict):
3180        alp = parmDict[hfx+'alpha']/refl[4+im]
3181        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
3182        return alp,bet
3183       
3184    hId = Histogram['hId']
3185    hfx = ':%d:'%(hId)
3186    bakType = calcControls[hfx+'bakType']
3187    fixedBkg = {i:Histogram['Background'][1][i] for i in Histogram['Background'][1] if i.startswith("_")} 
3188    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x,fixedBkg)
3189    yc = np.zeros_like(yb)
3190    cw = np.diff(ma.getdata(x))
3191    cw = np.append(cw,cw[-1])
3192       
3193    if 'C' in calcControls[hfx+'histType']:   
3194        shl = max(parmDict[hfx+'SH/L'],0.002)
3195        Ka2 = False
3196        if hfx+'Lam1' in (parmDict.keys()):
3197            wave = parmDict[hfx+'Lam1']
3198            Ka2 = True
3199            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3200            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3201        else:
3202            wave = parmDict[hfx+'Lam']
3203    else:
3204        shl = 0.
3205    for phase in Histogram['Reflection Lists']:
3206        refDict = Histogram['Reflection Lists'][phase]
3207        if phase not in Phases:     #skips deleted or renamed phases silently!
3208            continue
3209        Phase = Phases[phase]
3210        if histogram and not histogram in Phase['Histograms']:
3211            continue
3212        pId = Phase['pId']
3213        pfx = '%d::'%(pId)
3214        phfx = '%d:%d:'%(pId,hId)
3215        hfx = ':%d:'%(hId)
3216        SGData = Phase['General']['SGData']
3217        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3218        im = 0
3219        if Phase['General'].get('Modulated',False):
3220            SSGData = Phase['General']['SSGData']
3221            im = 1  #offset in SS reflection list
3222            #??
3223        Dij = GetDij(phfx,SGData,parmDict)
3224        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]  #TODO: need to do someting if Dij << 0.
3225        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3226        if np.any(np.diag(G)<0.) or np.any(np.isnan(A)):
3227            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
3228        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3229        Vst = np.sqrt(nl.det(G))    #V*
3230        if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
3231            if im:
3232                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
3233            elif parmDict[pfx+'isMag'] and 'N' in calcControls[hfx+'histType']:
3234                MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3235            else:
3236                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3237        badPeak = False
3238        # test to see if we are using multiprocessing here
3239        useMP,ncores = G2mp.InitMP()
3240        if len(refDict['RefList']) < 100: useMP = False       
3241        if useMP: # multiprocessing: create a set of initialized Python processes
3242            MPpool = mp.Pool(ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
3243            profArgs = [[] for i in range(ncores)]
3244        if 'C' in calcControls[hfx+'histType']:
3245            for iref,refl in enumerate(refDict['RefList']):
3246                if im:
3247                    h,k,l,m = refl[:4]
3248                else:
3249                    h,k,l = refl[:3]
3250                Uniq = np.inner(refl[:3],SGMT)
3251                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
3252                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
3253                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3254                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3255                refl[11+im] *= Vst*Lorenz
3256                 
3257                if Phase['General'].get('doPawley'):
3258                    try:
3259                        if im:
3260                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3261                        else:
3262                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3263                        refl[9+im] = parmDict[pInd]
3264                    except KeyError:
3265#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3266                        continue
3267                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
3268                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3269                iFin = np.searchsorted(x,refl[5+im]+fmax)
3270                if not iBeg+iFin:       #peak below low limit - skip peak
3271                    continue
3272                elif not iBeg-iFin:     #peak above high limit - done
3273                    break
3274                elif iBeg > iFin:   #bad peak coeff - skip
3275                    badPeak = True
3276                    continue
3277                if useMP:
3278                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin,1.))
3279                else:
3280                    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
3281                if Ka2:
3282                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3283                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
3284                    iBeg = np.searchsorted(x,pos2-fmin)
3285                    iFin = np.searchsorted(x,pos2+fmax)
3286                    if not iBeg+iFin:       #peak below low limit - skip peak
3287                        continue
3288                    elif not iBeg-iFin:     #peak above high limit - done
3289                        return yc,yb
3290                    elif iBeg > iFin:   #bad peak coeff - skip
3291                        continue
3292                    if useMP:
3293                        profArgs[iref%ncores].append((pos2,refl,iBeg,iFin,kRatio))
3294                    else:
3295                        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
3296        elif 'T' in calcControls[hfx+'histType']:
3297            for iref,refl in enumerate(refDict['RefList']):
3298                if im:
3299                    h,k,l,m = refl[:4]
3300                else:
3301                    h,k,l = refl[:3]
3302                Uniq = np.inner(refl[:3],SGMT)
3303                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position - #TODO - what about tabluated offset?
3304                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
3305#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
3306                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3307                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)             #TODO - skip if alp, bet tabulated?
3308                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3309                refl[11+im] *= Vst*Lorenz
3310                if Phase['General'].get('doPawley'):
3311                    try:
3312                        if im:
3313                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3314                        else:
3315                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3316                        refl[9+im] = parmDict[pInd]
3317                    except KeyError:
3318#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3319                        continue
3320                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
3321                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3322                iFin = np.searchsorted(x,refl[5+im]+fmax)
3323                if not iBeg+iFin:       #peak below low limit - skip peak
3324                    continue
3325                elif not iBeg-iFin:     #peak above high limit - done
3326                    break
3327                elif iBeg > iFin:   #bad peak coeff - skip
3328                    badPeak = True
3329                    continue
3330                if useMP:
3331                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin))
3332                else:
3333                    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]
3334#        print 'profile calc time: %.3fs'%(time.time()-time0)
3335        if useMP and 'C' in calcControls[hfx+'histType']:
3336            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfCW,profArgs):
3337                yc += y
3338            MPpool.terminate()
3339        elif useMP:
3340            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfTOF,profArgs):
3341                yc += y
334