source: trunk/GSASIIstrMath.py @ 3870

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

still mag incomm. wrong

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