source: trunk/GSASIIstrMath.py @ 3861

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

another try at incomm. mag. str. fctr.

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