source: trunk/GSASIIstrMath.py @ 3801

Last change on this file since 3801 was 3801, checked in by vondreele, 4 years ago

combine cell2AB & Gmat2AB in G2lat
fixes to mag. str. fctr. math. Correct sym. op. transformation of moments & fix moment to cartesian axis problem

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