source: trunk/GSASIIstrMath.py @ 4995

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