source: trunk/GSASIIstrMath.py @ 4997

Last change on this file since 4997 was 4997, checked in by vondreele, 3 months ago

move imports to top of G2ctrlGUi code so debugger can find it
fix FillUnitCell? bug when incommensurate mag.
latest revision of incomm. mag. str. factor stuff

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