source: trunk/GSASIIstrMath.py @ 4488

Last change on this file since 4488 was 4488, checked in by vondreele, 18 months ago

another try at incommensurate magnetic SF; fix bug in plotting mag moments.

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