source: trunk/GSASIIstrMath.py @ 4055

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

minor changes to incomm. mag. srt. fctr. stuff. Still wrong.

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