source: trunk/GSASIIstrMath.py @ 4185

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

mag math again

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