source: trunk/GSASIIstrMath.py @ 5015

Last change on this file since 5015 was 5015, checked in by vondreele, 2 years ago

changes to comments

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